UWO-TH-08/10 



OO 

o 

O 

(N 



Dynamics in the quantum Hall effect and the phase diagram of graphene 

E.V. Gorbar, 1 ^ V.P. Gusynin, 1 ^ V.A. Miransky, 2 -@ and LA. Shovkovy 3 '! 

1 Bogolyubov Institute for Theoretical Physics, 03680, Kiev, Ukraine 
2 Department of Applied Mathematics, University of Western Ontario, London, Ontario N6A 5B7, Canada 
3 Physics Department, Western Illinois University, Macomb, Illinois 61455, USA 

(Dated: August 28, 2008) 

The dynamics responsible for lifting the degeneracy of the Landau levels in the quantum Hall 
(QH) effect in graphene is studied by utilizing a low-energy effective model with a contact interac- 
tion. A detailed analysis of the solutions of the gap equation for Dirac quasiparticles is performed 
at both zero and nonzero temperatures. The characteristic feature of the solutions is that the order 
parameters connected with the QH ferromagnetism and magnetic catalysis scenarios necessarily co- 
exist. The solutions reproduce correctly the experimentally observed novel QH plateaus in graphene 
^ ' in strong magnetic fields. The phase diagram of this system in the plane of temperature and elec- 

tron chemical potential is analyzed. The phase transitions corresponding to the transitions between 
different QH plateaus in graphene are described. 
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I. INTRODUCTION 

In this paper, we analyze the dynamics in quantum Hall (QH) effect in graphene, a single atomic layer of graphite^ 
As was experimentally discovered in Refs. and theoretically predicted in Refs. 011111, an anomalous quantization 
takes place in this case: the filling factors are v = ±4(n + 1/2), where n = 0, 1, 2, ... is the Landau level index. For 
each QH state, a four-fold (spin and sublattice-valley) degeneracy takes place. These properties of the QH effect are 
intimately connected with relativistic-like features in the graphene dynamics^ ^ 10 ' 11 ^ 12 

In recent experiments ; 13 i 14 it has been observed that in a strong enough magnetic field, B > 20 T, the new QH 
plateaus with v — 0, ±1 and ±4 occur. This is attributed to the magnetic field induced splitting of the n = and 
ri = 1 Landau levels (LLs). It is noticeable that while the degeneracy of the lowest LL (LLL), n = 0, is completely 
lifted, only the spin degeneracy of the n = 1 LL is removed. 
£3 ', O n theoretical side, there are now two leading scenarios for the description of these plateaus. One of them is the 
l—— '■ QH ferromagnetism (QHF) i 15 i 16 ' 17 i 18 ' 19 (The dynamics of a Zeeman spin splitting enhancement considered in Ref. I20I 
j is intimately connected with the QHF.) The second one is the magnetic catalysis (MC) scenario in which Dirac 
^ . masses are spontaneously produced as a result of the excitonic condensationi 21 i 22 i 23 i 24 For a brief review of these two 
\Q ' scenarios, see Ref. di| 

, The QHF scenario is connected with the theory of exchange-driven spin splitting of Landau levels 2 ^ and utilizes 
OO ■ the dynamical framework developed for bilayer QH systems^ The underlying physics relies on the fact that the spin 
| and/or valley degeneracy of the one-particle states is lifted by the repulsive Coulomb interaction in a many-body 
\Q ■ system at half filling. The argument is the same as that behind the Hund's rule in atomic physics. The Coulomb 
energy of the system is lowered by antisymmetrizing the coordinate part of the many-body wave function. Because of 
OO the Fermi statistics of the charge carriers, the corresponding lowest energy state must be symmetric in the spin-valley 

' degrees of freedom, i.e., it is spin and/or valley polarized. 
J> , On the other hand, the MC scenario is based on the phenomenon of an enhancement of the density of states 
in infrared by a strong magnetic field, which catalyzes electron-hole pairing (leading to excitonic condensates) in 
relativistic-like systems. The essence of the MC phenomenon is the dimensional reduction D — > D — 2 in the pairing 
dynamics on the LLL with energy E = (containing both electron and hole states). In two dimensions, this reduction 
implies a non-zero, proportional to \eB\/2irhc, density of states in infrared. The latter is responsible for a Cooper- 
like electron- hole pairing even at the weakest attractive interaction between electrons and holes. This universal 
phenomenon was revealed in Ref. [28| and was first considered in graphite in Refs. I^fiol. 

The difference between the QHF and MC scenarios is in utilizing different order parameters in breaking an ap- 
proximate J7(4) symmetry of the Hamiltonian of graphene. This symmetry operates in the sublattice-valley and spin 
spaces. While the QHF order parameters are described by densities of the conserved charges connected with diagonal 
generators of the non-Abelian subgroup SU(4) C C7(4), the order parameters in the MC scenario are Dirac mass 
terms. 

One may think that the QHF and MC order parameters should compete with each other. However, as was recently 
pointed out by three of the authors^ the situation is quite different: these two sets of the order parameters necessarily 
coexist, which implies that they have the same dynamical origin. The physics underlying their coexistence is specific 
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for relativistic-like dynamics that makes the QH dynamics of the U (4) breakdown in graphene to be quite different 
from that in the bilayer QH systems 2 ! whose dynamics has no relativistic-like features. 

The main goal of this paper is a detailed study of the dynamics responsible for lifting the degeneracy of the Landau 
levels in the quantum Hall effect in graphene using the model of Ref. |2^. To get the benchmark results that are 
unobscured by the various types of possible disorder ) 30 ' 31 i 32 i 33 the analysis in this study is done for graphene in the 
clean limit. By taking into account a considerable improvement in samples quality seen in graphene suspended above 
a Si/SiC>2 gate electrode 3 -^ or above a graphite substrate^ it is expected that the clean limit already provides a 
reasonable qualitative description for some real devices (the role of disorder in this dynamics will be briefly considered 
in Sec. ED) 

The main tool in our analysis is a gap equation for the propagator of Dirac quasiparticles. The highlights of the 
analysis are as follows: 

1. The coexistence of the QHF and MC order parameters is a robust phenomenon, which is mostly based on the 
kinematic and symmetric properties of the QH dynamics in graphene. 

2. The process of filling the LLs is described by varying the electron chemical potential fiQ. A set of the solutions 
of the gap equation at a fixed fio is quite rich. The stable solution is selected as the solution with the lowest 
free energy density. The obtained results for the QH effect qualitatively agree with the experimental data in 
Refs. 0QJ. 

3. The existence of two types of the Dirac masses in the QH dynamics in graphene is established. Both of them 
play an important role in the dynamics. 

4. The phase diagram in the plane of temperature T and electron chemical potential /iq is analyzed. The phase 
transitions corresponding to the transitions between different QH plateaus are described. 

The paper is organized as follows. In Sec. |TT] we start by describing the general features of the model itself as well 
as the many-body approximation used in its analysis. After that, in Sec. 11111 we derive the gap equation for Dirac 
quasiparticles in graphene at zero and nonzero temperatures. The necessity of the coexistence of the QHF and MC 
order parameters in the solutions of the gap equation is shown. The analysis of the quasiparticle dynamics at the 
LLL, which is relevant to the v = 0, ±1 QH plateaus, is presented in Sec. HVl There we first give a detailed derivation 
of the analytic results of Ref. at zero temperature. Then, we consider the nonzero temperature case by utilizing 
numerical calculations. In a similar fashion, in Sec. |V1 the quasiparticle dynamics at the n = 1 Landau level is 
analyzed. In Sec. IVI1 we summarize our findings in the form of the phase diagram of graphene in the T — (Uq plane. 
The obtained phase diagram is rich and it allows to describe all the recently discovered novel plateaus (as well as 
the plateaus v = ±3 and v = ±5 which have not been observed yet) in graphene in strong magnetic fields . 13 i 14 We 
also discuss the correspondence between our results and the experimental data and point out that the coexistence 
of the QHF and MC order parameters could have important consequences for edge states, whose relevance for the 
dynamics in graphene has been recently discussed in Refs. [20ll36ll37l . Detailed derivations of some key results used in 
our analysis are presented in four Appendices at the end of the paper. 

II. MODEL: GENERAL DESCRIPTION 

A. Model: Hamiltonian and gap equation 

Our approach is based on the gap equation for the propagator of Dirac quasiparticles. For the description of the 
dynamics in graphene, we will use the model introduced recently in Ref. |23 . which in turn is a modification of the 
model in Refs. I9l ll0ll2ll . Let us start from the description of the latter. In this model, while quasiparticles are confined 
to a 2-dimensional plane, the electromagnetic (Coulomb) interaction between them is three-dimensional in nature. 
The low-energy quasiparticles excitations in graphene are conveniently described in terms of a four-component Dirac 
spinor = (tpKAs, 4>KBs, ipK'Bs, tpK'As) which combines the Bloch states with spin indices s = ± on the two different 
sublatticcs (A, B) of the hexagonal graphene lattice and with momenta near the two inequivalent valley points (K, 
K') of the two-dimensional Brillouin zone. The free quasiparticle Hamiltonian can be recast in a relativistic-like form 
with the Fermi velocity vf ~ 10 6 m/s playing the role of the speed of light: 

H = v F J d 2 r*( 7 1 7 r x + 7 2 7Ty) (1) 

where r = (x, y) is the position vector in the plane of graphene and ^ = vfrty is the Dirac conjugated spinor. In Eq. (JTJ) , 
j u with v = 0, 1, 2 are 4x4 gamma matrices belonging to a reducible representation of the Dirac algebra, namely, 
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FIG. 1: The diagrammatic form of the gap equation in the Hartree-Fock (mean field) approximation. The upper (lower) 
diagram corresponds to the form of the gap equation with the long-range Coulomb (contact) interaction. The indices denote 
quasiparticle spins. 



7" = f 3 <£> (t 3 , ir 2 , —ir 1 ), where the Pauli matrices f 4 and t % , with i = 1,2, 3, act in the subspaces of the valleys (K, 
K') and sublattices (A, B), respectively^ The matrices satisfy the usual anticommutation relations {7^,7"} = 2g tlv , 
where g^ v = diag (1, — 1, —1) and (1, v = 0, 1, 2. The canonical momentum 7t = (ir x ,ir y ) = —ifr'V + eA/c includes the 
vector potential A corresponding to a magnetic field Bj_ , which is the component of the external magnetic field B 
orthogonal to the ccy-plane of graphene. 

The Coulomb interaction term has the form 

H c = ^ J d 2 rd 2 r'^(r)-$(r)Uc(r-r')¥(r')-$(r'), (2) 

where E/c( r ) is the Coulomb potential in a magnetic field, calculated in the random phase approximation (RPA) in 
Ref. [13, see Eq. (46) there. The Hamiltonian H = Hq + He possesses a global U{4) symmetry discussed in the next 
subsection. The electron chemical potential //o is introduced by adding the term — /io^v]/ to the Hamiltonian density. 
This term also preserves the C7 (4) symmetry. The Zeeman interaction is included by adding the term hbB^ct 3 ^, 
where B = |B| and fiB = eTi/(2mc) is the Bohr magneton. Here we took into account that the Lande factor for 
graphene is — 2 (our convention is e > 0). The spin matrix a 3 has eigenvalue +1 (—1) for the states with the 
spin directed along (against) the magnetic field BiS Such states will be called spin up (down) states. Because of the 
Zeeman term, the Z7 (4) symmetry is broken down to a symmetry U(2) + x C/(2)_, where the subscript ± labels the 
spin of the states on which this subgroup operates (see the next subsection). 

The dynamics will be treated in the Hartree-Fock (mean field) approximation, which is conventional and appropriate 
in this case^ i 10 ' 15 i 16 i 21 Then, at zero temperature and in the clean limit (no impurities), the gap equation takes the 
form: 

G-^u, u') = S-^u, u') + ihj°G{u, u')~f°5(t - t')U c (r - r') - th^tr [y°G(u, u)] S 3 (u - u')U { c F) {0), (3) 

where u = (t,r), t is the time coordinate, U c (0) is the Fourier transform of Uc(r) at k = 0, G{u,u') = 
?i~ 1 (0|T\E , (u)^r('ti')|0) is the full quasiparticle propagator, and 

iS-^^u') = [(ihdt + - ti B B<7 3 )-f - v F (n ■ j)] S 3 (u-u') (4) 

is the inverse bare quasiparticle propagator. Note that while the second term on the right hand side of Eq. J3]) describes 
the exchange interaction, the third one is the Hartree term describing the direct interaction. The diagrammatic form 
of the gap equation is shown in Fig. [lja). 

As will be shown in Sec. IIIII below, in order to determine all the order parameters, the analysis of the gap equation 
([3]) has to be done beyond the LLL approximation, which is a formidable problem. Because of this, we follow the 
approach of Ref. I29I and replace the Coulomb potential J7c( r ) in the gap equation by the contact interaction Gj nt (5 2 (r). 
Thus, we arrive at 

G^(u, «') = 5 _1 (tt, u') + ihG in a°G(u, u)-f°6 3 (u - u') - ihG in a° tr[ 7 °G(u, u)]6 3 (u - u'), (5) 

where G- ln t is a dimensionful coupling constant. As we will see later, in the analysis it would be more convenient to use 
a dimensionless coupling constant A = Gi nt A/(4:7T 3 ^ 2 h 2 Vp) instead of Gj nt (note that A is the energy cutoff parameter 
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which is required when a contact interaction is used) . The corresponding diagrammatic form of the equation is given 
in Fig.[TJb). A similar approximation is commonly used in Quantum Chromodynamics (QCD), where the long range 
gluon interaction is replaced by the contact Nambu-Jona-Lasinio one. This leads to a good description of many 
features of the nonperturbative dynamics in low energy region of QCD (for a review see, for example, Ref. l40h . By 
taking into account the universality of the MC phenomenon and the fact that the symmetry and kinematic structures 
of equations j3| and ((5]) are the same, we expect that approximate gap equation j5]) should be at least qualitatively 
reliable for the description of the LLL and the first few LLs, say, with n = ±1. 

As to the value of the cutoff A, note that, in a strong magnetic field, the Landau scale 

e B ee ^2fi|eB ± |u!,/c~ 424V|£i[T]| K (6) 

is the only relevant energy scale in the dynamics with the Coulomb interaction. This suggests that the ultraviolet 
cutoff A should be taken of order cb in the approximation with the contact interaction. The dimensionful coupling 
constant Gi n t then becomes Gi nt ~ 4ir 3 ^ 2 h 2 VpX/eB- 

Before concluding this section, the following remark concerning the present approximation is in order. While there 
is Debye screening at nonzero chemical potential, the situation is more complicated near the Dirac point with /i = 0. 
In that case, while for subcritical values of the Coulomb coupling constant^ the polarization effects lead only to its 
screening without changing the form of the Coulomb interactions at large distances^ they lead to a drastic change 
of the form of the interactions for a supercritical coupling^ 3 - In the present work, the dynamics with a subcritical 
coupling is utilized, when no dynamical gaps are generated without a magnetic field (this is in agreement with the 
experiments^). In our approximation, utilizing smeared contact interactions with an ultraviolet cutoff A ~ es, the 
contribution of large energies u> > £b is suppressed much stronger than for the subcritical Coulomb like interactions. 
However, because the dominant contribution in the gap equation comes from energies u> <C es, we expect that the 
present approximation is qualitatively reliable even near the Dirac point. 



B. Model: Symmetry and order parameters 

The Hamiltonian H = Ho + He, with Ho and He given in Eqs. (fTl ) and $Z§, respectively, possesses the U(4) 
symmetry with the following 16 generators (see for example Refs. [il)ll2l"l ): 

<7 a a a a a o~ a 

-7T®k, 7F®^ 3 ' ^T®7 5 , and — ® 7 3 7 5 , (7) 
I Li I A 

where I4 is the 4x4 Dirac unit matrix and a a , with a = 0, 1, 2, 3, are four Pauli matrices connected with the spin 
degrees of freedom (cr° is the 2x2 unit matrix). In the representation used in the present paper (for the definition 
of the 7" matrices, see the previous subsection), the Dirac matrices 7 3 and 7 s = i7°7 1 7 2 7 3 arc 

where I is the 2x2 unit matrix. Note that while the Dirac matrices 7 and 7 = (7 1 ,7 2 ) anticommute with 7 s and 
7 5 , they commute with the diagonal matrix 7 3 7 5 = — 7 5 7 3 , 

7V=(J_° 7 ). (9) 

The matrix 7 3 7' 5 is called a pseudospin operator. 
The total Hamiltonian 

fltat =H + J d 2 r [/i B S*V 3 * - /i * f *] (10) 

possesses a lower symmetry. Because of the Zeeman term iibBWcf^, the E/(4) symmetry is broken down to the 
"flavor" symmetry U{2) + x C/(2)_, where the subscript ± corresponds to spin up and spin down states, respectively. 
The generators of the U(2) s , with s = ±, are h ® P s , —ij 3 ® P s , 7 s <£> P s , and 7 3 7 5 ® P s , where P± = (1 ± cr 3 )/2 are 
the projectors on spin up and down states. 

Our goal is to find all solutions of Eq. ([5]) both with intact and spontaneously broken SU (2) s symmetry, where 
SU(2) S is the largest non-Abelian subgroup of the U(2) s . The Dirac mass term As^Ps^ = A s ^-f Ps^, where 
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A s is a Dirac gap (mass)j^ is assigned to the triplet representation of the SU(2) S , and the generation of such a 
mass would lead to a spontaneous breakdown of the flavor SU(2) S symmetry down to the U(l) s with the generator 
7 3 7 5 <g> P ai & 10 i 21 There is also a Dirac mass term of the form A S ^7 3 7 5 P S ^]> that is a singlet with respect to SU(2) S , 
and therefore its generation would not break this symmetry. On the other hand, while the triplet mass term is even 
under time reversal T, the singlet mass term is T-odd (for a recent review of the transformation properties of different 
mass terms in graphene, see Ref. EH ) . Note that the possibility of a singlet Dirac mass like A was first discussed in 
relation to graphite about 20 years ago.— 

The masses A s and A s are related to the MC order parameters (^ , 7 3 7 5 P S ^I') and (^fP s ^f). In terms of the Bloch 
components of the spinors, the corresponding operators take the following forms: 

A s : *7 3 7 5 P S * = ^ KAs 1pKAs ~ ^ K , As ^K'As ~ ^kb^KBs + ^ K , Bs ^K'Bs, (11) 

A s : *P S * = 1pl <As 1pKAs + 1pK, As 1pK>As ~ i'KBs^KBs - Bs i'K'B s - (12) 

The expressions on the right hand side further clarify the physical meaning of the Dirac mass parameters as the 
Lagrange multipliers that control various density imbalances of electrons at the two valleys and the two sublattices. 
In particular, the order parameter (|12[) . connected with the triplet Dirac mass, describes the charge density imbalance 
between the two sublattices, i.e., a charge density wave^2i 

As revealed in Ref. and will be discussed in detail in the next section, these MC order parameters necessarily 
coexist with QHF ones in the solutions of the gap equation ((5]). More precisely, for a fixed spin, the full inverse 
quasiparticlc propagator takes the following general form [compare with Eq. Q]: 

iGj\u,u') = [(ihdt + us + As7 3 7 5 )7° -uf(tt-7)- A s + A s7 3 7 5 ] <5 3 (u - u'), (13) 

where the parameters fi s , fl s , A s , and A s are determined from gap equation ([5]). Note that the full electron chemical 
potentials \i± include the Zeeman energy =pZ with 

Z ~ I-J-bB = 0.67P[T] K. (14) 

The chemical potential jx s is related to the density of the conserved pseudospin charge ^j 3 ^ P s ^, which is assigned 
to the triplet representation of the SU(2) S . Therefore, unlike the masses A s and A s , the chemical potentials ^3 = 
(/.t+ — /i-)/2 and p, s are related to the conventional QHF order parameters: the spin density (vf , t<7 3 \I r ) and the 
pseudospin density (^E , ^7 3 7 5 P S \I'), respectively. In terms of the Bloch components, the corresponding operators take 
the following forms: 

*v* = i E (v'L+vw-C-Vva-), (15) 

" k=K,K' a=A,B 

% \ * t 7 3 7 5 P s * = i> ] KAs ^KAs ~ ^ K , As 1pK'As + ^KBsi'KBs - ^ K , Bs ^K'Bs- (16) 

By comparing the last expression with Eq. (fl~2"l) , we see that while the triplet MC order parameter related to A s 
describes the charge density imbalance between the two graphene sublattices, the pseudospin density (related to p, s ) 
describes the charge density imbalance between the two valley points in the Brillouin zone. On the other hand, as 
seen from Eq. (|15p . ^,3 is related to the conventional ferromagnetic order parameter (\E''cr 3 ^) . 

The following remark is in order. Because of the relation 7 s = i7°7 1 7 2 7 3 , the operator in Eq. (| 16[) can be 
rewritten as i\E , 7 1 7 2 P s \I/. The latter has the same form as the anomalous magnetic moment operator in Quantum 
Electrodynamics (QED). However, unlike QED, in graphene, it describes not the polarization of the spin degrees of 
freedom but that of the pseudospin ones, related to the valleys and sublattices. Because of that, this operator can be 
called the anomalous magnetic pseudomoment operator. 

Let us describe the breakdown of the £7(4) symmetry down to the U{2) + x U(2)_ flavor symmetry, responsible for 
a spin gap, in more detail. Because of the Zeeman term, this breakdown is not spontaneous but explicit. The point 
however is that, as was shown in Ref. I20I, a magnetic field leads to a strong enhancement of the spin gap in graphene. 
Such an enhancement is reflected in a large chemical potential /13 = (fi + — /i-)/2 3> Z and the corresponding QHF 
order parameter (^V 3 ^). But as was pointed out in Ref. [2^ and will be shown below in Sec. IIV1 it is not all. There 
is also a large contribution to the spin gap connected with the flavor singlet Dirac mass A3 = (A + — A_)/2 and the 
corresponding MC order parameter (^ , 7 3 7 5 cr 3 \E'). This feature leads to important consequences for the dynamics of 
edge states in graphene (see Sees. HVIand IVI[). 

As will be shown in Subscc. flV C[ the spin gap may remain large even in the limit when the Zeeman energy Z = /i^B 
goes to zero. In this limit, a genuine spontaneous breakdown of the f/(4) takes place. In the realistic case with a 
nonzero but small Z, one can say that a quasi-spontaneous breakdown of the U (4) is realized. 
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The 17(2)+ x f/(2)_ is an exact symmetry of the total Hamiltonian H tot (fT0|) of the continuum effective theory. 
However, as was pointed out in Ref. [l?] (see also Refs. [l9|, [22|, and l46h. it is not exact for the Hamiltonian on the 
graphenc lattice. In fact there are small on-site repulsion interaction terms which break the £7(2)+ x £/(2)_ symmetry 
down to a U(l)+ x Z-2+ X £7(1)_ x Z2- subgroup, where the elements of the discrete group Z^s are 7 s (g) P a + 14 ® 
and the unit matrix. Unlike a spontaneous breakdown of continuous symmetries, a spontaneous breakdown of the 
discrete symmetry Z2±, with the order parameters (tyP±'ty) and (\t r '''7 3 7 5 .P±\E f ), is not forbidden by the Mermin- 
Wagner theorem at finite temperatures in a planar system.— This observation is of relevance for the description of 
the ground state responsible for the v = ±1 plateaus (see Subscc. HVD[) . 

Thus, there are six order parameters describing the breakdown of the £7(4) symmetry: the two singlet order 
parameters connected with fi 3 and A 3 and the four triplet ones connected with fx± and A±. 

By extracting the location of the poles in full propagator G(u, u'), which is given in terms of the sum over separate LL 
contributions in Eq. (|A27[) in Appendix [XJ it is straightforward to derive the dispersion relations for the quasiparticles 
in graphene. The dispersion relations for LLs with n > 1 are 

w<# = -li a + ajis ± yf ne 2 B + (A s + aA s f , (17) 

where a = ±1 and the two signs in front of the square root correspond to the energy levels above and below the Dirac 
point. In the case of the LLL, which is special, the corresponding dispersion relations read 

ljW = -fx s + a U s sign(eB ± ) + A.) + A, sign(eB ± ). (18) 

As shown in Subsec. IA 21 in Appendix \X[ the parameter a in Eqs. (jTTJ) and ([T5)l is connected with the eigenvalues of 
the diagonal pseudospin matrix 7375 in Eq. ©. For the LLs with n > 1, the value a = ±1 in (jTTJ) corresponds to 
the eigenvalues =pl of 7 3 7 5 . On the other hand, for LLL, the value a = ±1 in (|18[) corresponds to sign(e_B+) x (=pl) , 
with =pl being the eigenvalues of 7 3 7 5 . 

One can see from Eqs. (|17|l and l|18p that at a fixed spin, the terms with a are responsible for splitting of LLs. We 
will return to this issue in Sec. IIVI 



III. GAP EQUATION: EXPLICIT FORM AT T = AND T AND COEXISTENCE OF QHF AND MC 

ORDER PARAMETERS 

In this section, we will present the explicit equations for the Dirac masses and the chemical potentials at zero and 
finite temperature. In particular, it will be shown that the QHF and MC order parameters necessarily coexist. 

The equations for the Dirac masses A s and A s and the chemical potentials ^ s and jl s follow from the matrix form 
of the gap equation in Eq. (fSJ) and expression (|13[) . Their derivation, while straightforward, is rather tedious. It is 
considered in Appendix lAl in detail. At zero temperature, the equations are 

A [ 

A s = — < - [sign(/z s - (j, s )9(\fj, s -fi s \- Eos) - si S n (Ms + Vs)0{\fJ- s + A* s | - E 0s )] sign(e J B_ L ) 



E 

n=0 



(A s + A s )8(E+-\ f i s -fi s \) , (A S -A S )6(E- S -\fi s + fi s \) 



^ n s 



E n 



l + 0(n-l) 



(19) 



A s = j < - [sign(/x s - fr s )9(\fj, s -fi s \- E£ s ) + sign(/i s + p, s )9(\fj, a + fi a \ - E 0s )] sign(e J B_ L ) 



E 



(A s + A S )6(E+ - \fi a - fi s \) (A s ~ A S )6(E- S - | Ms + fi„\) 



E + 



Ens 



[l + e(n-i)}}, 



(20) 



.4 



(A S + A S )6(E+ -\fi a -fi a \) (A, - A a )6(Eo 8 - \fi s + fi a 



E 0s E 0s 



sign(e£_L) 



J2 [-sign(/i a - fi a )0(\n. -fi a \- E+ s ) + sign( Ms + (M s )9(\(x a + fi„\ - E~)] [1 + 6(n (21) 



7 




(A. + A S )9(E+ 



jS a |) {A B -A a )9{Eo„-\n. + fi.\) 



Os 

E 0s E 0s 



sign(eB_L) 



+ E [agnO*. - - -S+J + sign^+^^dMs+^l --B ns )] [l + 0(n-l)] L 



(22) 



n=0 



where the step function is defined by 9(x) = 1 for x > and 0(x) = for x < 0. Regarding the other notation, 
fi± = fio =F Z is the bare electron chemical potential which includes the Zeeman energy Z = fisB, and E^ s = 

\J ne 2 B + (A~ s ± A s ) 2 are quasiparticle energies. In these equations, we introduced a new energy scale, A, that plays 
an important role throughout the analysis. It is determined by the value of the magnetic field and the coupling 
constant strength, 



.4 



8ttHc 4A 

The second term on the right hand side in Eq. ([22]) is defined as follows: 

x = E x « 



(23) 



(24) 



where 



-2Al 



(A S +A S )6(E+ 



\H. - fx s \) (A, - A s )0(££- S - | Ms + 



E, 



n.s 



sign(eBj.) 



+ E [ si S n (^ " - A.I - + sign(/i s + A.)^(lA*- + h\ - Ks)] [1 + *(n - 1)] 



(25) 



n=0 



The following comment is in order here. Because of the Hartree term in the gap equation the equations for the 
spin up and spin down parameters do not decouple: they are mixed via the X term in Eq. (|22[) . Fortunately, it is 
the only place affected by the Hartree term. As shown in Appendix [B] this fact strongly simplifies the analysis of 
the system of equations (fT9")) - (|2"2"|) . This point also clearly reflects the essential difference between the roles played 
by the exchange and Hartree interactions in the quasiparticle dynamics of graphene. While the former dominates 
in producing the QHF and MC order parameters, the latter participates only in the renormalization of the electron 
chemical potential, which is relevant for the filling of LLs. 

Since the step functions in the above set of equations depend on [a s ± jl s and A s ± A s , it is more convenient to 
rewrite the gap equations for the following set of parameters 



Ai±) = A s ± A s 



(±) - u ±u 



(26) 



In the numerical analysis, we always consider a nonzero temperature. This is implemented by utilizing the Matsubara 
formalism. Using the identities 



T E 

n— — oc 

CO 

T E 



1 



smh(E/T) 



[{2n + l)irT + i^i} 2 + E 2 

-i(2n+ 1)ttT + fx 
[(2n + 1)ttT + ifj] 2 + E 2 



2E cosh{E/T) + cosh(/x/T) ' 



sinh(/i/T) 



2cosh(£;/T) + cosh(/i/T)' 



(27) 
(28) 



it is straightforward to write the equations at nonzero temperature. One can check that the prescription for modifying 
Eqs. (fl9 ]) -([22] ) at T ^ is to replace 



S ign(^)9(\^\-EZ 
6{E 



J ns 



l/4 T) l) 



sinh ; 



,(±) 



cosh -^p + cosh 



sinh 



cosh -jj^ + cosh ^jr- 



(29) 
(30) 
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This leads to the following set of equations: 



(±) 



/'■s 



A/ 2 (AW, M ( ± ))+2A/ 2 (A( ± ), / 
where A^ and fig are given in Eq. (|26[) . and 



,(T) 



2A/ 2 fA^,^ ) 



(31) 
(32) 



h (aw 



,(T) 



sinh f A y ^ 


— sj_ sinh (t 5- ) 


cosh ^-^f^- 


) + cosh 



E 



2Ai ±} sinh (^p) 




-Ens 


cosh ^-^p 


+ cosh 







s j_ sinh ( 



A 



(±) 



sinh 



cosh 



cosh 



n=l COSh 



2 sinh ( ^jr- 



cosh 



(33) 



(34) 



with s± = sign(e_B^) and E. 



i 

ns 



Ai ±: 



Let us now show that the QHF and MC order parameters should always coexist in this dynamics. Suppose that 
Eqs. (|3"Tj) and (f3"2"| have a solution with some of the chemical potentials fj,f being nonzero but the Dirac masses 
being zero, A s = 0. Then, the left hand side of Eq. (|31[) is equal to zero. On the other hand, taking into account 
expression (f3"3")) for the function /i, we find that for As = the right hand side of this equation takes the form 



/i(o,/4 T) ) = 



-s i sinh 



1 + cosh 



,(=F) 



= — s I tanh 



2T 



(35) 



and it could be zero only if all chemical potentials /x s disappear, in contradiction with our assumption. Therefore 
we conclude that the QHF and MC order parameters in this dynamics necessarily coexist indeed. This is perhaps one 
of the central observations in this study. 

Which factors underlie this feature of the graphene dynamics in a magnetic field? It is the relativistic nature of the 
free Hamiltonian Hq in Eq. ([1]) and the special features of the LLs associated with it. To see this, note that while the 
triplet Dirac mass A s multiplies the unit Dirac matrix I4, the triplet chemical potential jl s comes with the matrix 
7 3 7 5 7° in the inverse propagator Gj 1 in Eq. (fl~3|) . Let us trace how these two structures are connected with each 
other. The point is that there are terms with i7 1 7 2 sign(ei?^) matrix in the expansion of the propagator G 8 over LLs 
[see Eq. (|A21[) in Appendix [A] . Taking into account the definition 7 s = i7°7 1 7 2 7 3 , we have i7 1 7 2 = 7 3 7 5 7°. Then, 
through the exchange term ~ 7°G S 7° in gap equation ([5]), the A s term in the inverse propagator Gj 1 necessarily 
induces the term with the chemical potential jl s . In the same way, the singlet Dirac mass A s in Gj 1 is connected 
with the singlet chemical potential fi s . 

These arguments are based on the kinematic structure of gap equation ([5]) , which is the same as that for equation 
([3]) with the Coulomb interaction. Taking into account the universality of the MC phenomenon, we conclude that the 
coexistence of the QHF and MC order parameters is a robust feature of the QH dynamics in graphene. 

The necessity of the coexistence of the QHF and MC order parameters can be clearly seen in the case of the 
dynamics on the LLL. As follows from Eq. (|A27|) in Appendix [5] the LLL propagator contains only the combinations 
— [i s + A s sign(e-Bj_) and /x s sign(e£?jj + A s . Therefore, in this case, the QHF and MC parameters not only coexist 
but they are not independent, which in turn reflects the fact that the sublattice and valley degrees of freedom are not 
independent on the LLL. In particular, by using Eqs. (fTTj) , (fl2"l) . (|15[) . and one can easily check that, because of 
the projector V- = [1 — i7 1 7 2 sign(eB_L)]/2 in the LLL propagator [see Eqs. (|A21[) and (|A22[) ]. the operators ^P s ^ 
and 4'7 3 7 5 P s v]> (ty^ 3 ^ 5 P s fy and ^PsVP), determining the order parameters related to (i s and A s (jj, s and A s ), coincide 
up to a sign factor sign(ei?jjj^i This fact in particular implies that in order to determine all the order parameters, it 
is necessary to analyze the gap equation beyond the LLL approximation. 

The important point, however, is that this special feature of the LLL takes place only on an infinite plane. In 
real graphene samples with boundaries the situation is different: the QHF and MC parameters on the LLL become 
independent. 49 ' 50 As is discussed in Sec. IVI1 this leads to important consequences for the dynamics of edge states on 
the LLL. 
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IV. DYNAMICS ON LLL: v = 0, v = ±1, AND v = ±2 PLATEAUS 

As was already discussed in Introduction, at magnetic fields B < 10 T, the plateaus with the filling factors 
1/ = ±4(n + 1/2) are observed in the QH effect in graphene^ At stronger magnetic fields, new plateaus, with v = 
and v = ±1 occur: while the former arises at B > 10 T, the latter appear at B > 20 T^i4 In this section, we will 
describe the dynamics underlying these new plateaus, and the plateaus v = ±2 corresponding to the gap between 
the LLL and the n = 1 LL, by using the solutions of the gap equation presented in the next subsection. We will 
consider positive v and fig (the dynamics with negative v and /io is related by electron-hole symmetry and will not 
be discussed separately) . As will be shown below, there is a large number of the solutions corresponding to the same 
fiQ. In order to find the most stable of them, we compare the free energy density Q for the solutions. The derivation 
of the expression for fl is presented in Appendix [Cl 



A. Overview of analytic solutions at LLL 



The v = 0, v = ±1 and v = ±2 plateaus are connected with a process of doping of the LLL, which is described 
by varying the electron chemical potential /io. Therefore we start our analysis by reviewing the solutions to the gap 
equations in the case when /io is much less than the Landau energy scale, i.e., /io < eb. At zero temperature the 
corresponding gap equations are analyzed analytically in Appendix [Bj It is concluded there that only the following 
three stable solutions arc realized: 

(i) The solution with singlet Dirac masses for both spin up and spin down quasiparticles, 

A+ = /i+ = 0, fx+ = fi + -A, A+ = s ± M, 

(36) 

A_ = p,_ = 0, fx- = p,- + A, A_ = -s ± M. 

[By definition M = A/(l — A) and A = 4AA/(y / 7reg), see Eq. (|B9[) and its derivation in Appendix [Bl] This solution is 
energetically most favorable for < /io < 2A+ZM- It is one of several solutions with nonvanishing singlet Dirac masses 
and we call it the 5*1 solution (here S stands for singlet). Because of the opposite signs of both the masses A + and 
A_ and the chemical potentials /i+ and /i_, the explicit breakdown of the f (4) symmetry down to U(2)+ x t/(2)_ 
by the Zeeman term is strongly enhanced by the dynamics. Since all triplet order parameters vanish, the flavor 
U{2) + x f/(2)_ symmetry is intact in the state described by this solution. As discussed in Subscc. flV CI below, the 
51 solution corresponds to the v — plateau. 

(ii) The hybrid solution with a triplet Dirac mass for spin up and a singlet Dirac mass for spin down quasiparticles, 

A+ = M, fi + =As ± , M +=M+-4A A + = 0, 
A_ =0, fi_ = 0, fj,_ = p,_ - 3A, A_ = -s ± M. 

It is most favorable for 2 A + Z < /io < 6 A + Z. Wc call it the HI solution (here H stands for hybrid, meaning that 
the solution is a mixture of the singlet and triplet parameters). In this case, while the SU(2)+ C U(2)+ symmetry 
connected with spin up is spontaneously broken down to U(l)+ (whose generator is j 3 j 5 (g)P + ), the SU(2)_ C U(2)_ 
symmetry connected with spin down remains intact. As will be shown in Subsec. IIVD1 the HI solution corresponds 
to the v = 1 plateau. 

(iii) The solution with equal singlet Dirac masses for both spin up and spin down quasiparticles 

A+ = /i+ = 0, n+ = /i+ - 7 A, A + = -s±M, 
A_ = /i_ =0, m_ = /2_ - 7 A, A_ = -s ± M. 

It is most favorable for /io > 6^4 + Z. We call it the S2 solution. (Note that the dynamics in the n = 1 LL will set an 
upper limit for the range where the S2 solution is the ground state, see Sec IVl below.) In the state given by the S2 
solution, the U(4) symmetry is broken down to U(2) + x t/(2)_ only by the Zeeman term. Indeed, the singlet masses 
and the dynamical contributions to the chemical potentials arc of the same sign for both spin orientations and thus 
have no effect on breaking any symmetry. As will be shown in Subsec. HVEl the S2 solution corresponds to the v = 2 
plateau connected with the gap between the filled LLL and the empty n = 1 LL. 

The free energy densities for the above three solutions are given by the following expressions (see Subsec. IB 61 in 
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FIG. 2: Free energy density versus the electron chemical potential /io for several different solutions, found analytically (left 
panel) and numerically (right panel) in a range of /io relevant to the dynamics in the lowest Landau level. The numerical results 
are shown for a nonzero but small temperature, T = 1 K. The values of the electron chemical potential are given in units of the 
Landau energy scale es, and the free energy densities are given in units of es/? 2 , where I = ^Hc/\eB± | is the magnetic length. 



Appendix [B| : 



n 
n 



2irhc 
\eB±\ 
2ttHc 
\eB±\ 
2nhc 



(M + A + 2Z + K), for < [i <2A + Z, 

(M - A + Z + h + n ), for 2A + Z < /i < 6A + Z, 

(M - 7A + h + 2/i ) , for 6A + Z<^ Q , 



(39) 
(40) 
(41) 



where the parameter h is defined in Eq. (|B56[) . We note that although the parameters of the solutions jump abruptly 
at the transition points, [io = 2A + Z and /io = GA + Z , their free energy densities match exactly. We conclude, 
therefore, that first order phase transitions take place at these values of the electron chemical potential /io. 

The free energy densities in Eqs. (|39p ~ (|41[) are shown as functions of the chemical potential ^iq in the left panel in 
Fig. [21 In order to plot the results, we took M — 4.84 x 10~ 2 es and A = 3.90 x 10~ 2 es which coincide with the 
values of the corresponding dynamical parameters in the numerical analysis. For comparison, the numerical results 
at nonzero but sufficiently small temperature are shown in the right panel of Fig. [21 As we see, the agreement is 
very good. It is interesting to note that the singlet-type numerical solution, given by the solid line, spans both the 
SI and S2 solutions, as well as the intermediate (metastable) branch connecting them. In addition to the 5*1, HI, 
and S2 solutions, numerical results for several other (metastable) solutions are shown. The metastable solutions are 
discussed in Subsec. TlV Fl below. 



B. Numerical analysis at LLL 



In this subsection, we give the key details regarding our numerical analysis. 

Throughout this paper the default choice of the magnetic field in the numerical calculations is B = 35 T. The 
corresponding Landau energy scale is £b|b=35 t ~ 2510 K. In order to do the numerical calculations in the model at 
hand, we use a simple rcgularization method that renders the formally defined divergent sum in Eq. (j3"3")l finite. In 
particular, we redefine the corresponding function as follows: 



A( 



sinh 



(- 



Y~) — s± sinh 



°° 2A.i ±) sinh(^) n(y^e B ,A) 



cosh (^rj + cosn 
where k{x, A) is a smooth cutoff function defined by 

sinh (A/<5A) 



n— 1 Ens 



cosh ( 



cosh 



k(x, A) 



cosh (x/5A) + cosh (A/5A) 



(42) 



(43) 
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with A = 5000 K and SA = A/20 = 250 K. The value of A corresponds to an approximate point of the high-energy 
cut-off, and the value of i5A gives the extent of the smearing region in either direction from A. (Note that the energy 
scale A is about the same as the energy of the n — 4 Landau level at B = 35 T.) 

One should emphasize that the specific choice of the cutoff energy scale A has little effect on the qualitative as well 
as quantitative results of our analysis, provided the dynamical energy scales A and M = A/(l — A) are kept fixed (see 
the discussion in the end of this subsection). Here we assume that the value of the cutoff is sufficiently large to avoid 
the reduction of the phase space relevant for the quasiparticlc dynamics at the n — and n = 1 LLs. 

Because of the cutoff function k(x, A) the sum over n on the right hand side of Eq. (|42|) is rapidly convergent. In the 
numerical calculations, therefore, a sufficiently good accuracy may be achieved by keeping a finite number of terms 
in the sum. The optimum choice for the maximum value of index n is n max = [l4A 2 /e^] , where the square brackets 
mean the integer number nearest to the result in the brackets. This choice is large enough to insure a high precision 
and, at the same time, it is small enough to make the calculation fast. 

While the /2-function in Eq. (|34[) is finite, for consistency we redefine it in the same way as function f\ by smoothly 
cutting off the contributions of large-n LLs, 



5± sinh(41)-sinh(^) _ - 2sinh(^).(V^,A) 



cosh (^^f^J + cosn \T^J n=1 cos h ( T~ > + cos ^ 

where k(x,A) is defined in Eq. (|43p . The numerical result for the sum in fa is also approximated by dropping the 
terms with n > n max where n max is given above. 

By analyzing the solutions to Eqs. (|3"T]) and (|32|) at very low temperatures, we reproduce all the analytic solutions 
derived in Appendix [BJ For the choice of the magnetic field B = 35 T the values of the two dynamical energy 
parameters A and M are given by 

A ps 98 K, M ps 122 K. (45) 

As is easy to check, these correspond to the dimensionless coupling A ps 0.196. Here one should keep in mind that the 
smooth-cutoff rcgularization used in our numerical calculations is not the same as in the analytical calculations [see, 
for example, Eq. (|B6[) in Appcndix[B]] Despite this difference, all analytical results agree very well even quantitatively 
with the corresponding numerical ones when expressed in terms of A and AI parameters. 



C. Plateau v = 



The plateau v = is connected with a range of electron chemical potentials in the vicinity of the Dirac neutral 
point with /zo = 0. In this case the 51 solution with singlet Dirac masses of opposite sign for spin up and spin 
down quasiparticles, see Eq. (|31))) . is most favorable energetically and therefore is the ground state solution, provided 
fio < 2 A + Z (other solutions related to the Dirac neutral point are discussed in Subsec. IIVFI below). 

From dispersion relation (|18[) . we find that while u;+ = —jiQ + Z + M + A is positive for spin up states, cj_ = 
— Ho — Z — AI — A is negative for spin down states, i.e., the LLL is half filled (the energy spectrum in this solution is 
a independent). Therefore there is a nonzero spin gap AEq = uj + — w_ associated with the v — plateau. The value 
of this gap is AE = 2(Z + A) + 2M. 

While no exact symmetry is broken in the state described by the 51 solution, the explicit spin symmetry breaking 
by the Zeeman term Z is strongly enhanced by the dynamical contribution M + A. In this case, it is appropriate to 
talk about the dynamical symmetry breaking of the approximate spin symmetry. This is also evident from studying 
the temperature dependence of the MC and QHF order parameters in Fig. [3J In the two panels, we compare the 
results in the models with the exact (left panel) and approximate (right panel) spin symmetry. In the first case we take 
Z = and see that the spontaneous spin-symmetry breaking occurs at low temperatures. The symmetry is restored 
at about T ps 0.9M in a typical second order phase transition (recall that we work in the mean-field approximation). 
In the second case, a nonzero Zeeman energy term (Z ps 23.51 K at B = 35 T) breaks the spin symmetry explicitly 
and its restoration is impossible even at very high temperatures. However, even in this latter case, there is a well 
pronounced crossover (around T ps 0.9M) between the regimes of low and high temperatures, which can be quantified 
by the relative strength of the bare Zeeman and dynamical contributions. 

The order parameters for the solution 51 versus the electron chemical potential /io are shown in Fig. 0] for several 
different values of the temperature. At T = this solution is the ground state for fj,o < 0.09es. At sufficiently low 
temperature, the main qualitative feature of this solution is that the singlet Dirac masses for spin-up and spin-down 
quasiparticles have opposite signs, A + = — A_. This defines the configuration of the MC order parameters that is 
formally invariant under the time reversal symmetry. (Of course, the time reversal symmetry is still explicitly broken 
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T/M T/M 

FIG. 3: Temperature dependence of the nontrivial order parameters in the v = QH state, described by the 5*1 solution. The 
results in a model with a vanishing Zeeman energy (Z — 0) are shown in the left panel, and the results in a realistic model 
with a nonzero Zeeman energy (Z ^ 0) are shown in the right panel. Note that p,± = A± = in both cases. The values of the 
temperature and the order parameters are given in units of the dynamical scale M. 
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FIG. 4: Order parameters for the singlet solution S1/S2 as 
values of temperature. 
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ions of the electron chemical potential for several different 



by the external magnetic field.) As the temperature increases, the approximate relation A + ss — A_ may hold at 
Ho w 0, but deviations from such a relation grow with increasing fi . 

It should be emphasized that the solution SI is continuously connected with the solution S2 responsible for the 
v = 2 QH plateau, see Subsec. TlV El below. At low temperatures, the intermediate branch between the SI and S2 
solutions is metastable. At high temperatures, however, it becomes stable and the qualitative difference between the 
two solutions disappears. 

The conclusion that the v = state is related to the spin gap agrees with the scenario in Ref. [20| and the experiments 
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reported in Refs. [l3l36l . The fact established in the present paper that both /X3 and the singlet Dirac mass A3 contribute 
to the gap AEq is noticeable. As was already pointed out in Sec. IIIH unlike the case of an infinite plane, in graphcnc 
samples with boundaries, the parameters /X3 and A3 are independent on the LLL. As will be discussed in Sec. I VII this 
fact could have important consequences for the dynamics of edge states. 

In conclusion, the following comment is in order. As one can see in the right panel in Fig[2j besides the 51 
solution, there is another, triplet (T), solution around the Dirac neutral point. In the T solution, given in Eq. (|B33p 
in Appendix [B] both spin up and spin down quasiparticle states have a triplet Dirac mass. Calculating the difference 
of the free energy densities for these two solutions, one finds that 8fl = ilsi — = —Z\eB\/nhc. Therefore, it is the 
Zeeman term which makes the SI solution more favorable: without it, the SI and T solutions would correspond to 
two degenerate ground states. It would be interesting to figure out the role of the small on-site repulsion interaction 
term a 17 i 19 i 22 ' 46 mentioned in Subsec. Ill Bl in choosing the genuine ground state in the present dynamics. 

D. Plateau v = 1 

As was pointed out in Subsec. IIV Al for larger /xq the hybrid HI solution (f3"T|) . with a triplet Dirac mass for 
spin up quasiparticles and a singlet Dirac mass for spin down quasiparticles, is most favorable. It is the ground 
state for 2 A + Z < fxo < 6 A + Z. As one can easily check by using Eq. (fP8|) . while now w+ > 0, the energies 
lu^_ ' and = uj^ ^ are negative. Consequently, the LLL is now three-quarter filled and, therefore, the gap 

A.Ei = wV - ' — lo^ ' = 2{M + A) corresponds to the v = 1 plateau. Notably, the Zeeman term docs not enter the 
value of the gap. Unlike the v — state, therefore, the gap in the v = 1 state is directly related to the spontaneous 
breakdown of the flavor symmetry 577(2)+. 

The last point regarding the nature of the ground state described by the HI solution has important consequences 
for the physical properties of the v = 1 QH state. Since the coupling constant Gi n t in the present model is proportional 
to 1/cb (see Subsec. Ill Ap . Eq. (f23|) implies that the dynamical parameters A and M, and therefore the gap AE\, 
scale with the magnetic field as y/\eB±\. This fact agrees with the dependence of the activation energy in the v = 1 
state observed in Ref. [bH . 

The critical temperature at which the 5C/(2)+ symmetry is restored, i.e., when the triplet parameters /i+ and A+ 
vanish, is T c ~ 0.9M ~ 110A'. The restoration is described by a conventional second order phase transition. 

The temperature dependence of the hybrid HI solution is rather interesting too. This is summarized in Fig. [5] 
where the nontrivial order parameters and chemical potentials are shown for several values of the temperature in the 
range from 1 K to 100 K. One of the most spectacular features of this dependence is a revival of the singlet mass A+ 
at finite temperature shown in Fig. [5] (recall that it vanishes at zero temperature). This phenomenon is intimately 
connected with the general conclusion in Sec. IIIII that at a fixed value of spin s and any value of temperature, there 
are no nontrivial solutions of the gap equation with the both masses A s and A s being zero. Indeed, at T > T c , when 
the triplet mass A + vanishes, the absence of the A + would contradict this conclusion (note that as Fig. [5] shows, the 
revival of this mass occurs even at subcritical T). Note also that in the case of spin down quasiparticles, the triplet 
parameters /i_ and A_ are identically zero but the singlet mass A_ remains nonzero at all temperatures. 

These results arc obtained in the mean field approximation and for the Hamiltonian iTtot (jlOP - which is symmetric 
under the U(2) + x [/(2)_. However, as was already pointed out in Sec. Ill Bl above, this symmetry is not exact 
for the Hamiltonian on the graphene lattice. In that case, it is replaced by U(l)+ x Z2+ x J7(l)_ x Z2-, where 
the elements of the discrete group Z<i± are 7 s ® P± + I4 ® Pip and the unit matrix. It is important that unlike a 
spontaneous breakdown of continuous symmetries, a spontaneous breakdown of the discrete symmetry Zi±, with the 
order parameters (fyP±fy) and ^f 3 j 5 P± 1 i>) , is not forbidden by the Mermin- Wagner theorem at finite temperatures 
in a planar system4i This point strongly suggests that there exists a genuine phase transition in temperature related 
to the v = 1 state in graphene. 

E. Plateau v = 2 

At zero temperature, the S2 solution (|38p with equal singlet Dirac masses for spin up and spin down states is most 
favorable for /iq > 6 A + Z. It is easy to check from Eq. (fT5)) that both w + and w_ are negative in this case, i.e., 
the LLL is completely filled. This solution corresponds to the v = 2 plateau when the value of the electron chemical 
potential is in the gap between the LLL and the n = 1 LL. 

The nonzero temperature results for the order parameters of the solution 5*2 versus the electron chemical potential 
/io are shown in Fig. [U At T = this solution is the ground state for fiQ > 0.24es. As we see, even at high 
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FIG. 5: Order parameters for the solution HI as functions of the electron chemical potential fio for several different values of 
temperature. 



temperatures, the MC order parameters satisfy the same approximate relation. A+ w A_. Such a configuration 
breaks neither spin nor sublattice-valley symmetry of graphene. 



F. Metastable solutions on LLL 



As was already pointed above, in addition to the three stable solutions 51, HI, and 52, describing the v = 0, 
v = ±1, and v = ±2 QH plateaus, the numerical analysis of the gap equations reveals other, metastable, solutions. 

One of such solutions is the T solution with nonzero triplet Dirac masses for both spin up and spin down quasipar- 
ticles. In the model of graphene used in this paper, the explicit analytical form of this solution is given in Eq. (|B33[) 
in Appendix [Bj Note that because there is a contribution of the bare Zeeman term Z oc eB in the gap AEi for this 
solution, the corresponding activation energy in the v = 1 state scales with eB differently from the ^/|ei?| law in the 
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FIG. 6: Nontrivial order parameters of metastable solutions H2 (left panel) and 53 (right panel) as functions of the electron 
chemical potential /io- In calculation, the temperature is taken nonzero but small, T — IK. The values of the electron chemical 
potential are given in units of the Landau energy scale es, while the order parameters are given in units of the dynamical scale 



hybrid HI solution. 

In addition to the triplet solution, there exist also metastable hybrid (H2) and singlet (53) solutions. Their free 
energy densities are shown in Fig. [5] together with the energy densities of the other solutions. As seen, neither the 
H2 solution nor the 53 one can have sufficiently low free energy density to become the genuine ground state. 

The following remark is in order. Unlike all the other solutions, the solutions H2 and 53 cannot be found analytically 
at T = 0, see Appendix [B] By making use of the numerical analysis, we find that these two extra solutions are such 
that /4 ~ ±E^ S . At exactly zero temperature, it is problematic to get such solutions analytically because Eqs. (fT§|) - 
([^)) contain undetermined values of the step functions, e.g., 9(1/1^ \—E^ s ). In contrast, at a nonzero temperature, the 

step functions are replaced by smooth expressions, see Eqs. (|29p and (|30p. and numerical solutions with fig ps ±E^ 
are easily found. The order parameters for the solutions H2 and 53 are shown in Fig. [6] 



V. DYNAMICS ON n = 1 LL 



In the previous section, we analyzed solutions of the gap equations under the condition that only states on the LLL 
can be filled, \fi s ± jl s \ <C £s = \J 2ft\eB ±]v 2 F / c. Since all the dynamically generated parameters are much less than 
es, this condition implies that the bare chemical potential fi also has to satisfy a similar inequality, /io < tg. In this 
section, we extend the analysis by considering the dynamics with being of the order of the Landau scale es, i.e., 
we study the regime when quasiparticle states on the first Landau level, n = 1 LL, can be filled. 



A. Analytic solutions at T = 

We will start from the gap equations at zero temperature, which are given in Eqs. (fTl?)) -(|22 p in Sec. IIII1 In order 
to get their solutions for /j,q ~ es, we will follow the same steps of the analysis as in Appendix [B] for the LLL. The 
corresponding analysis for the n = 1 LL, including the calculation of the free energy density for the solutions, is done 
in Appendix [Dl It is shown there that the following five stable solutions are realized (see the end of Subsec. fD 3|) : 

(f-i) The singlet type solution (f-I-f-I) (here / stands for first; the nomenclature used for the n = 1 LL solutions is 
defined in Appendix [Dj : 

A+=A+=0, At+=A 2 + -7A, A+ = -s ± M, 

(46) 

A_ = /2_ = 0, /i_ = pj— - 7A, A_ = -sj_M 
coincides with the solution S2 given by Eq. (|38| in the analysis of the LLL. It takes place for QA + Z < fj,o < 
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7 A + \J ( 2 B + M 2 — Z, and its free energy density is 



Q = -^^(M + 2Mo-7A + fr), (47) 

where h is given in Eq. (|B56[) . According to Subsec. [TV El this solution corresponds to the regime with the filled 
LLL and the empty n = 1 LL and is connected with the v = 2 plateau. 

(f-ii) The hybrid type solution (f-I-f-II) 

A+ = m+ = 0, /Lt + = fa - 11A, A+ = -a x M, 

- M — Mi , , M + Mi (48) 

A_ = — -, = -Asx, M- = M- - 10A, A_ = -s_l — -, 



with Mi given in Eq. (|D5]l in Appendix[Ql takes place for 9A + yJt 2 B + Mj 2 - Z < ^ < H-A + \As + A/2 ~ z > 
and its free energy density is 

\eBx\ ( 3M + Mi ir „ „ 3fc + fci\ 

where /ii is given in Eq. (|E)22|) . As is shown in Subsec . IV Bl below, this solution corresponds to the v = 3 plateau, 
(f-iii) The singlet type solution (f-I-f-III) 



A+ = fa = 0, /i+ = fa - 15A, A + = -sx M, 
A_ = fa =0, //_=//_- 13A, A_ = -sx Mi 



(50) 



is realized for 13^4 + y/e B ~~+ M 2 — Z < fiQ < 15 A + \] t 2 B + M 2 + and its free energy density is 



leBxl /M + Mi , h + h!\ 

= -V^ T" — + 4 ^o - 27 A - 2e B + 2Z + — — ^ . (51) 



As is discussed in Subsec. IV Bl this solution corresponds to the v = 4 plateau. 

(f-iv) The hybrid type solution (f-II-f-III) 

- M-Mi . a a M + Mi 

A+ = , fx + = -Asx, M+ = M+ - 18A, A + = -s_l , 

A_ = /2_ = 0, [i- = fa - 17 A, A_ = -sx Mi 



takes place for 17A + \/e 2 B + M 2 + Z < (j,q < 19A + \Jt 2 B + M 2 + Z, and its free energy density 



is 



\eBx\ fZM 1 + M A 3h! + h\ 
= -V^ ~ A + 5 t i -A3A-3e B + Z+^- . (53) 



2nhc V 

This solution corresponds to the v = 5 plateau (see Subsec. fVBp . 
(f-v) The singlet type solution (f-III-f-III) 

A+ = fa = 0, =fa- 21A, A+ = -sx Mi 



(54) 

A-=fa = 0, /i_ = fa- 21A, A_ = -s_l Mi 



is realized for jiQ > 21A + \/e B + M 2 + and its free energy density is 



fi = -J|^(Mi + 6/i -63A-4e s + /i 1 ). (55) 
27mc 
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FIG. 7: Nontrivial order parameters of the S-type numerical solution that contains the analytical solutions (f-i), (f-iii) and (f-v) 
as parts, connected by two intermediate solutions. 



This solution corresponds to the v = 6 plateau connected with the gap between the filled n = 1 LL and the empty 
n = 2 LL. 

It should be emphasized that the above analytical solutions do not cover the whole range of the values of the 
electron chemical potential around the n = 1 LL. In particular, there are no analytical solutions found in the following 
four intervals: 



7A 
llA-r 
15A4 
19,4 H 



M 2 - Z < fj, < 9 A 



M 2 - Z, 



M 2 -Z<n < 13A 



Ml - Z, 



M 2 



Z < /i < 17 A 



M\ + Z, 



M 2 + Z < fiQ < 21 A + Je 2 B + M 2 + Z 



(56) 
(57) 
(58) 
(59) 



The difficulty in finding analytical solutions at T = on these intervals is related to the ambiguities in the definition 
of some step functions in gap equations (TT9"]) - (|2"2"1) . The same problem, albeit in a weaker form, was also encountered 
in the analysis of dynamics at the LLL (see Subscc. [lVF[) . As in that case, we remove the ambiguities by considering 
a nonzero temperature case. The results at T = can then be obtained by taking the limit T — > 0. The details of 
our numerical analysis are given in the next subsection. 



B. Numerical analysis, n — 1 LL 



By performing a nonzero temperature analysis numerically, wc find that the solutions (f-i), (f-iii), and (f-v), found 
analytically, are in fact continuously connected. They are parts of a more general solution S (here S stands for singlet) 
that exists at all values of /io- At small and intermediate values of [iq, this solution includes solutions 51 and 52, see 
Fig. |U At larger values of fig, relevant for the dynamics of n = 1 LL, the solution S is shown in Fig. [71 

As seen in Fig. [JJ the solution S consists of five pieces defined on five adjacent intervals of (j,q. Three of them are 
the analytical solutions (f-i), (f-iii), and (f-v), as defined in the previous subsection. Their intervals of existence are 
Mo/ £ s ~ 1-27, 1.5 < fio/^B ^5 1-6 and Ho/tB ^ 1.83, respectively. These intervals are in agreement with the analytical 
results if one takes Mi r* 111 K, or in terms of the Landau energy scale, Mi = 4.42 x 10~ 2 eb- The other two pieces 
of the solution S extend the singlct-typc analytical solution to the intermediate intervals. 

At T = the solution S describes the ground state in exactly the same regions of validity that are found analytically 
for solutions (f-i), (f-iii), and (f-v) in the previous subsection. This can be concluded from the energy consideration: 
among all numerical solutions the parts of the solution S have the lowest free energy density there. Analyzing the 
quasiparticle spectra by using the dispersion relation in Eq. (|17[) . we find that the solutions (f-i), (f-iii), and (f-v) 
describe the v = 2, v = 4, and v = 6 QH states, respectively. 

From the symmetry viewpoint, none of the three parts of the singlet solution break any exact symmetries in the 
model. However, the part (f-iii) of the solution, describing the v — 4 QH state, corresponds to a quasi-spontaneous 
breakdown of the ?7(4) symmetry down to the U(2)+ x J7(2)_. Indeed, by using Eq. (jTTJ) , one can check that the 
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FIG. 8: Nontrivial order parameters of the extended hybrid solution (f-ii) which determines the ground state for the v = 3 QH 
plateau in graphene. 



LLL is half filled and the energy gap between the pairs of the pseudospin degenerate spin-up and spin-down states of 
the n = 1 LL is given by AE^ ~ 2(Z + A) + {M 2 — M\)j2e B . As we see, the spin splitting by the Zeeman term 2Z 
is strongly enhanced by the dynamical contribution 2A. 

This is somewhat similar to the enhancement of the spin splitting in the v = QH state, discussed in Subsec. IIV CI 
However, there is an important qualitative difference between the cases of the LLL and the n = 1 LL: It is only the 
dynamical contribution to the chemical potentials (but not the Dirac masses) that substantially affects the splitting in 
the v = 4 QH state. Indeed, the dynamical contribution due to the Dirac masses in the gap AE4, i.e., (A/ 2 — Mf)/2eB, 
is very small because M ~ Mi <C e B -) As a result, the gap AE4 is substantially smaller than the LLL spin gap AEq 
(AE a < AE /2). 

Because of having nonvanishing triplet order parameters in the extended hybrid solutions (f-ii) and (f-iv) , the flavor 
U{2) + x U(2)_ symmetry of graphene is partially broken in the corresponding ground states. By using dispersion 
relation (|17[) in the analysis of the quasiparticle spectra, we find that these solutions describe the v = 3 and v = 5 
plateaus corresponding to the quarter and three-quarter filled n = 1 LL, respectively. In the case of the extended 
solution (f-ii), the spin-down flavor subgroup SU(2)_ C U(2)- is broken down to [/(1)_, while the spin-up flavor 
subgroup U(2)+ is intact. Similarly, in the case of the extended solution (f-iv), the spin- up flavor subgroup SU(2) + C 
U{2) + is broken down to U(l)+, while the spin-down flavor subgroup E/(2)_ is intact. Up to small corrections due 
nonzero Dirac masses, the energy gaps AE3 and AE5 associated with the (f-ii) and (f-iv) solutions are equal to 2A. 
Note that these gaps are substantially smaller than the LLL gap AE\ (AE3, AE5 < AE\/2). 

The analytical hybrid solutions (f-ii) and (f-iv) get continuous extensions to the left and to the right from their regions 
of validity found analytically in the previous subsection. In fact, they extend all the way to cover the neighboring 
"forbidden" regions defined in Eqs. (|55|) (15^1) . The first two "forbidden" interval are covered by the extension of the 
solution (f-ii) to the interval 7 'A + \J + M 2 — Z < fj,o < 13A + \J f 2 B + M\ — Z. The non-trivial Dirac masses 
and chemical potentials for this numerical solution are shown in Fig. [8l The last two "forbidden" intervals, see 
Eqs. ([58]) and (|59|) . are covered by the extension of the solution (f-iv) to the interval 15^4 + ^/e 2 B + AP + Z < fi < 
21A + \/e B + M 2 + Z. The non-trivial parameters for this solution are shown in Fig. [5] 

In fact, the extended solutions (f-ii) and (f-iv) are the ground states in their whole regions of existence. This is seen 
in Fig. [ini where we plot the difference between the free energy density of the hybrid type solutions and the singlet 
one. The results for the extended hybrid solutions (f-ii) and (f-iv) are shown by the solid line and the long-dashed 
line, respectively. 

In Fig. [10] we also show the results for another hybrid solution that was found numerically. It exists in the interval 
of no that could potentially be relevant for the v = 4 QH state. However, its free energy density is higher than that 
for the solution S 7 and therefore it is unstable. 

With increasing the temperature, we find that the extended hybrid solutions (f-ii) and (f-iv) responsible for the 
v = 3 and v = 5 QH states gradually vanish. Their regions of existence shrink and their free energy densities approach 
the free energy density of the singlet solution S. At temperatures above T„~^ ~ ri^ -5 ' ~ 0.4M ~ /2, they 

cease to exist altogether, and the ground state is described by the singlet solution which does not break any exact 
symmetries of the model. 
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FIG. 9: Nontrivial order parameters of the extended hybrid solution (f-iv) which determines the ground state for the v = 5 QH 
plateau in graphene. 
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FigH 




VI. DISCUSSION: PHASE DIAGRAM, EXPERIMENT, DISORDER, AND EDGE STATES 

By summarizing the numerical results for the ground states at different temperatures, we obtain the phase diagram 
of graphene in the plane of temperature T and electron chemical potential (1q shown in Fig. 1111 The areas highlighted 
in green correspond to hybrid solutions with a lowered symmetry in the ground state. These regions are separated 
from the rest of the diagram by phase transitions. At the boundary of the u = 1 region, the transition is of first 
order at low temperatures and of second order at higher temperatures. The transitions to/from the QH v = 3 and 
v = 5 states are of second order. It should be kept in mind, however, that here the analysis is done in the mean-field 
approximation and in a model with a simplified contact interaction. Therefore, the predicted types of the phase 
transitions may not be reliable. In particular, the contributions of collective excitations, which are beyond the mean- 
field approximation, may change the transitions to first order type. Also, the types of the transitions may be affected 
by the inclusion of disorder and a more realistic long-range Coulomb interaction. Despite the model limitations, we 
still expect that Fig. [TT] correctly represents the key qualitative features of the phase diagram of graphene at least in 
the case of the highest quality samples. 

In Fig. [TT] the regions highlighted in blue correspond to the ground states with a quasi-spontaneous breakdown of 
the spin symmetry. In the case of the LLL and the n = 1 LL, such are the v = and v = 4 QH states, in which the 
quasi-spontaneous breakdown of the approximate U(4) symmetry down to J7(2)+ x C/(2)_ is enhanced by dynamical 
contributions. Because of the explicit breakdown by the Zeeman term, there is no well-defined order parameter 
associated with this symmetry breakdown. Also, there is no well-defined boundary of the corresponding regions in 
the diagram. In the plot, this feature is represented by the fading shades of blue at the approximate boundaries of 
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FIG. 11: Schematic phase diagram of graphene in the plane of temperature and electron chemical potential. The values of 
chemical potential are given in units of the Landau energy scale es, and the values of temperature are given in units of the 
dynamical scale M. 



the v = and v = 4 regions. 

As considered in detail in Sec. IIV1 the physical properties of the v = and v = 1 QH states arc determined by 
the dynamics of the LLL. The corresponding values of the gaps, AEo = 2(Z + A + M) and AE\ — 2(A + M), are 
largely determined by the dynamical contributions A and M of about equal magnitude. These two contributions are 
associated with the QHF and MC order parameters, respectively. 

The results of this study for the LLL at least qualitatively agree with the experimental data. 13 ' 14 By taking the 
dimensionlcss coupling A = 4AA/ (y/ire^) to be a free parameter and utilizing the cutoff A to be of the order of the 
Landau scale es, we arrive at the following scaling relations: A ~ \^J\eB±] and M ~ Xy/\eB±\. This implies the 
same type of scaling for the gap, AEi = 2(A + M) ~ Xy / \eB±\, associated with the u = ±1 plateaus. [Recently, the 
square root scaling of the activation energy in the v = 1 QH state was also obtained in the large- N approximation in 
Ref. HH.] By making use of our results, we find that the experimental value AEi ~ 100 K for B± = 30 T from Ref. [lil 
corresponds to A ~ 0.02. This estimate, however, should be taken with great caution: Because interactions with 
impurities are ignored and no disorder of any type is accounted for in the present model, it may not be unreasonable 
to assume that actual values of A are up to an order of magnitude larger. 

As to the n = 1 LL, we found that there are the gaps AE3 = AE$ ~ 2A and AE4 ~ 2(Z + A) corresponding to 
the plateaus v = 3, 5 and v = 4, respectively [the contributions of Dirac masses are suppressed by a factor of order 
(M/es) 2 there]. Therefore the gaps AE3 = AE5 and AE4 are mostly due to the QHF type order parameters and are 
about a factor of two smaller than the LLL gaps AEi and AEo, respectively. On the other hand, the experimental 
data yield AE4 ~ 2Z, and no gaps AE3, AE$ have been observed. 13 ' 14 We believe that a probable explanation of 
this discrepancy is that, unlike Z, the value of the dynamically generated parameter A corresponding to the |n| > 1 
LLs will be strongly reduced if a considerable broadening of higher LLs in a magnetic field is taken into account 3i If 
so, the gap AE4 will be reduced to 2Z, while the gaps AE3 and AE$ will become unobservable. 

In order to estimate the value of a magnetic field at which the plateaus v = 3 and 5 could become observable, one 
can use the following arguments. Recently, in Ref. HH, a large width Y\ of 400 K was determined for the n = 1 LL. On 
the other hand, the plateaus v = 3, 5 could become observable if the gaps AE3 = AE$ ~ 2A calculated in the clean 
limit are at least of order Ti or larger^ The LLL gap AE\ ~ 100 K at \B_\_\ = 30 T corresponds to AE 3:5 < 50 K. 
Then, taking a conservative estimate T\ = 100 K and using A ~ y^eBjJ, we conclude that to observe the v = 3, 5 
plateaus, the magnetic fields should be at least as large as B ~ 100 T. 

Here it is also appropriate to mention the dynamics of sap less edge states, whose importance for the physics of 
the v = plateau has been recently discussed in Refs. I20ll36ll37l . Generalizing the analysis in Ref. HH, it has been 
recently foun d 49 ' 50 that for the SI solution (|36|) with the zigzag boundary conditions, such states exist only when 
the full Zeeman energy — /U-)/2 = Z + A is larger than the Dirac mass A± = M (at an armchair edge, gapless 
edge states exist for any value of a singlet Dirac mass). Because of that, for A smaller than 1, we find from the 
constraint Z > \A/(l — A) in the solution (i) and Eq. (|23[) with A ~ that the gapless edge states exist when 
\B ± \ > B { ™ ] - 8 x 10 4 A 4 /(1 - A) 2 T. Then, for the values of the dimensionless coupling A in the range 0.02 < A < 0.2, 
we find that 0.01 T < B^ < 200 T. As we see, B± is very sensitive to the choice of A. Therefore, in order to 

(cr) 

fix the critical value B ± more accurately, one should first utilize more realistic models of graphene that incorporate 
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disorder among other thingsi 31 i 32 This is a topic for future studies however. 

These results are of interest in connection with the interpretation of the v = Hall plateau. Indeed, the gapless edge 
states should play an important role in transport properties of graphcnc in a strong magnetic field. Their presence 
is expected to make graphene a so-called quantum Hall metal, while their absence should make it an insulatori 20 i 36 
The actual temperature dependence of the longitudinal resistivity at the v = plateau in Refs. Il3ll36l is consistent 
with the metal type. This conclusion may be disputed in view of the recent data from Ref. [3?] that reveal a clear 
plateau at v — 0, but the temperature dependence of the diagonal component of the resistivity signals a crossover to 
an insulating state in high fields. The latter observations do not seem to support the existence of gapless edge states. 

The analysis in this paper as well as in Refs. EH3 

suggests that the conditions for the existence and absence of 
gapless edge states depend sensitively on the type of the boundary conditions and the values of QHF and MC order 
parameters that characterize the nature of the corresponding QH state. Therefore, the dynamics of the edge states is 
very likely to be rich and full of surprises. 

In conclusion, we have shown that the QHF and MC order parameters in graphene are two sides of the same coin 
and they necessarily coexist. This feature could have important consequences for the QH dynamics, in particular, for 
edge states. The present model leads to a reasonable and consistent description of the new QH plateaus in graphene in 
strong magnetic fields. It would be desirable to extend the present analysis to a more realistic model setup, including 
the Coulomb interaction between quasiparticles, the quasiparticlc width, and various types of disorder. 
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APPENDIX A: QUASIPARTICLE PROPAGATOR AND THE GAP EQUATION 
1. Quasiparticle propagator: Expansion over LLs 

In this Appendix, the units with h = 1 and c = 1 are used. 

The full propagator G s {u,u') that corresponds to the inverse propagator in Eq. (|13|) is given by the following 
expression: 



G s (u,u') = i(u\ (id t + ^s)T -v F (n -7) +i/2»7 7 + iA^V 7 - A * \ u ') 
= i(u\ (id t + i u s )7°-ViK^-7)+^ s 7 1 7 2 -«A s 7°7 1 7 2 +A s 
({idt + j u s )7° - v f (tt • 7) + i/i»7 1 7 2 + *'A,7V7 2 - £») 
v f (tt ■ 7) + ifisl 1 ! 2 — iA a 7 7 7 + A 



i(u\ (id t + Msb ~ v f (tt ■ 7) + «'As7 1 7 2 - «A S 7°7 1 7 2 + A, 



l«'> 



(id t + Ms) 2 - vW + 2ipt s (idt + Ms)7°7 1 7 2 + 2zA s A s7 7 1 7 2 



-ieB ± vWr + K - A 2 - A; 



(Al) 



where u — (t,r) and r = (x, y). Our aim is to get an expression for this propagator as an expansion over LLs. For 
the Fourier transform in time we need to calculate 



G s (lo- r, r') = i[W- v F (TV r ■ 7)] (r| (M - v 2 F ir 2 - ieB ± v 2 F 7 V) 1 |r' 



(A2) 
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where the matrices W and M. are 

W = (w + A ; s )7 + ^ s 7 1 7 2 -*A s 7 7 1 7 2 + A fl , (A3) 
M = (w + ft +iMV7 2 ) 2 -(A s -iA s 7°7V) 2 - (A4) 

The operator 7r 2 has well known eigenvalues (2n + l)|eBj_| with n = 0, 1, 2, . . . and its normalized wave functions in 
the Landau gauge A = (0,B±x) are 



where H n {x) are the Hermite polynomials and / = y^c/|eSjJ is the magnetic length. These wave functions satisfy 
the conditions of normalizability 



d 2 rijj* nv {v)il> n , v ,(r) = 6 nn ,6(p-p'), (A6) 
and completeness 

oo 

oo 

£ / « p (r)V„ P (r') = £(r-r'). (A7) 

"=°-oo 

Using the spectral expansion of the unit operator (|A7[) . we can write 

i mxa 22 • n 2 l 2\— l i n 1 f (r-r') 2 .(a; + a;')(y - y')\ 
(r\(M-v 2 F iv 2 -ieB ± v 2 F ~/W) IO = ^ cxp (- ^ ^ 

°° 1 f ( r ~ r ') 2 \ 

X ^ A4 - (2n + l)v 2 F \eB ± \ - %v 2 F eB^ 2 Ln \2P ) ' (A8) 

where we integrated over the quantum number p by making use of the formula 7.378 in Ref. l53l . 

e- x2 H m {x + y)H n {x + z)dx = 2 n 7: 1 / 2 m\z n - m L^ m (-2yz), (A9) 



assuming m < n. Here L" are the generalized Laguerre polynomials, and L n = The matrix iv F eB±j 1 j 2 has 
eigenvalues ±v F \eB±\ 1 and thus one can write 

MO P-MO + p+MO 



X - (2n + l)v||eBx| - w|e J B_ L7 1 7 2 X - (2n + l)u|.|eBx| + u||eSx| - (2n + l)v|.|eB±| - w||&B±| ' 

(A10) 

where the variable £ and the projectors P± are 

P± = \ [I±i7'7 2 sign(e-Bl)] ■ (A12) 
Now, by redefining n — > n — 1 in the second term in Eq. (|A10[) , equality (|A8|) can be rewritten as 

n=0 *" 1 

where L_x = by definition and the phase 



$(r,r>) = - {x + X '^ y,) =-e 
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appears because in the presence of a constant magnetic field, the commutative group of translations is replaced by the 
noncommutative group of magnetic translations^ (note that the integration in Eq. (|A14[) is taken along the straight 
line). This implies that it has a universal character. By noting that 



7r x e = e 



Tr^e* - e **^, + £ZJ^, (A16) 

we see that propagator (|A2|) can be presented in the form of a product of the phase factor and a translation invariant 
part G s (ui;r- r'), 

G,(u; r, r') = e^'^G.fo r - r'), (A17) 

where 



G s (uj;r - r') = i 



W - vp-y 1 ( -id x - V —f ) - v Fl z ( -id y + 



2l 2 r 1 V y 21 



e 



2ttZ 2 ^ M- 2nvl\eB x \ 



It is important to emphasize that the phase factor does not affect the gap equation ([5]) because the latter contains 
the full propagator only at u' = u. 

The Fourier transform of the translation invariant part of propagator (|A18[) can be evaluated by first performing 
the integration over the angle, 



2jv 

ikr cos 8 



d6e lkrcos0 = 2nJ (kr), (A19) 



where Jq(x) is the Bcssel function, and then using the formula 7.421.1 in Ref . l53l. 

»-»-"!. (l^) Jbto)* = ^e-A«"L„ (j^) , (A20) 
valid for y > and Re a > 0. The result is given by 

with 

Dn.Cw, k) = 2W [V-L n (2fc 2 / 2 ) - (2fc 2 / 2 )] + 4^(k ■ (2fc 2 Z 2 ) , L° x = 0, (A22) 

describing the nth Landau level contribution (compare with corresponding expression for the standard Dirac propa- 
gator in Ref. HI). 



2. Equations for Dirac masses and chemical potentials 

In order to derive Eqs. (fllH) - (f2"2")) for masses and chemical potentials, we need to know the full propagator at u' = u, 
G s (u,u')\ u=u > = G s (u,u). As follows from Eq. (|A18[) . it is 

oo oo 

°^ U)= J 2^'° ) = 2^£ J 2^ W M-2nvi\eB ± \ - (A23) 

In what follows, it is convenient to work with eigenvectors of the matrices j lr y 2 and 7 . Since (7 1 7 2 ) 2 = — 1, the 
eigenvectors |si2) of the matrix 7 1 7 2 correspond to imaginary eigenvalues isn = ±i, i.e., 

7 1 7 2 ki2) = isi 2 |si 2 ) • (A24) 
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Similarly, since (7 ) 2 = 1, the eigenvectors \sq) of the matrix 7 correspond to eigenvalues so = ±1, i.e., 

7°l*o> = *o|*o>. (A25) 

Because 7 and 7*7 2 commute, we can use states |si 2 so) which are simultaneously eigenvectors of 7 1 7 2 and 7 with 
eigenvalues isi 2 and So, respectively. The vectors |si 2 so) form a complete basis. Therefore, any 4x4 matrix O can 
be represented as 



°= E °<^ So \s' l A){^ \ 



(A26) 



Now, taking into account that propagator (|A23j) contains only the unit, 7 , 7 1 7 2 , and -y^ 1 ^ 2 matrices [see Eqs. (|A3j) 
and (j A4|) ] . its expansion in the form (|A26[) has only diagonal terms with s' 12 = si 2 and s' — sq. Therefore, we can 
rewrite it as follows: 



G s (u,u) 



4tt2 2 



S12.S0. 



2tt 



(iO + Us ~ VsS12Sq)s + A s + A s Si 2 5q 

^0 ( w + ^ s - Ms s i2 s o) 2 - (A s + A s si 2 s ) 2 - 2v 2 F \eBx\n 



x {1 + si 2 sign(eB ± ) + [1 - si 2 sign(e J B_ L )] 9(n - 1)} \s 12 s ) (s 12 s \ 



(A27) 



The zeros of the denominator in the integrand define the dispersion relations for the Landau levels. In the case of 
n > 1, they are given by 



fi s + ap, s ± J2v 2 F \eB ± \n+ (A s + aA s ), 



(A28) 



where a = si 2 so (i-e., a = ±1) and the two signs in front of the square root correspond to the energy levels above and 
below the Dirac point. The case of the LLL is special because the numerator in the n = term in Eq. (|A27|) coincides 
with one of the zeros in the denominator. After taking this into account, we find the following dispersion relation: 



-fi s + a\jx s sign(eB ± ) + A s ] + A s sign(e J B_ L ). 



(A29) 



Note that the parameter a — ±1 in Eqs. (|A28[) and (|A29[) is connected with the eigenvalues of the diagonal pscudospin 
matrix 7375 in Eq. ©. Indeed, from the expression 7 s = i7°7 1 7 2 7 3 , one gets 7 3 7 5 = ij°j lr y 2 , i.e., the eigenvalues of 
7 3 7 5 arc — sqs\2- It is now easy to check for higher LLs that a = ±1 in Eq. (|A28j) corresponds to the eigenvalues =f1 
of 7 3 7 5 . On the other hand, as follows from Eq. (|A21|) . si 2 = sign(e£>jj on the LLL, and we find that in this case 
(j = ±1 corresponds to sign(eS^) x (^fl), with =pl being the eigenvalues of 7 3 7 5 . 
Integrating over w in Eq. (|A27j) . we obtain 



G, 



si2,s n=0 \ 



s sign(/i s - fl s s 1 2So)6Qfi s - fi s s 12 s \ - E° s ) 

(A s + A a s 12 so)0(EZ s - \fi s - ^ S 5 12 s |) \ 
E° j 

^ns J 

x {l + s 12 sign(eSx) + [l-si2sign( eJ B_L)]e(n-l)} |s 12 s )(si2So|, (A30) 

where E°„ = ^2v 2 F \eB ± \n + (A s + crA s ) 2 . 

Using this expression and the inverse bare and full propagators in Eqs. (|4]) and (|13p . we arrive at the following form 
of gap equation ([5|): 

00 

-/jl s s + /x s s 12 + AsS 12 s Q + A s = - ft s s + [l + si 2 sign(eBj_) + (1 - si 2 sign(eSj_)) 0(n - 1) 

n=0 



-s sign(^ s - /2 s si 2 s o )0(|/i a - p, a s 12 s \ - E° 



(A s + As 12 s )e{E^ s - \iis - £ s si 2 5 |) 
E a 



As °E E E I 1 + ^MeB±) + (1 - s' 12 sign(eB ± )) 6(n - 1) 

(A s/ + A s .s' 12 s' )0(E° s , - |/v - /W 12 s l 



-s sign(^ s - - fx sf s' u s' )e(\fi s ' - p, s >s' 12 s' \ - E° B ,) + 



E a , 

ns 



(A31) 
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where A = G; n t/(87rZ 2 ). The last term on the right-hand side of Eq. (|A31[) proportional to so is the Hartree contribu- 
tion. Finally, multiplying (|A31[) by 1, S12S0J s i2 5 an d sg, respectively, and taking the sum over s 12 and Sq, we obtain 
Eqs. CSJ-J22II. 



APPENDIX B: ANALYTIC SOLUTIONS OF GAP EQUATION FOR LLL AT T = 



In order to solve Eqs. (p~9|) - (| 2 1 [) for A s , A s , fl s as functions of fj, s , note that these equations contain 6- functions 
whose arguments suggest that the following three cases have to be considered: 

1. l^sT/M < |A S ± A s |; 

2. \fx s - /2 S | > |A S + A s |, \fi a +fx s \< \A S - A s | or \fx a - p,„\ < \A„ + A s \, \/j, s + p, a \ > \A S - A s |; 

3. |M s T/i s | > |A S ±A S |. 



1. The first case 



For |/z s =F Ms I < |A S ± A s |, the gap equations for Dirac masses take the form 



A s + A s = (Bl) 



71—0 



n s 



A - A 

A s -A s = f - [1 + g(n - 1)] ■ (B2) 

,1=0 Ans 

Equations for A s + A s and A s — A s are equivalent and since each equation admits both positive and negative solutions 
with the same absolute value, we have 

A s + A s =±(A S - A s ). (B3) 
This implies that one of the following should be true 

(a) A s = 0, or (b) A s = . (B4) 
Then, the gap equation for the nonvanishing parameter A s (or A s ) takes the form 

A S =AJ2^= S . [l + e(n-l)\. (B5) 
n=o ne% + A 2 

Let us first consider the case (a) and show that Eq. (|B5p can be equivalently represented in the following integral 
form: 

A, = ^r^ [ . tl (| 9 ), (Be, 

where A is a high energy cut-off up to which the low-energy effective theory is valid. After taking into account the 
identity coth(e^y/2) = 1 + 2 J^Li e~ yntB , we can integrate over y in Eq. (|B6|) by using the following table integral: 

%L e -y(nel+Al) „ f°° dy_ e -v(n£> B +&*) = (B7) 

/a 2 Vv Jo Vv 



where we replaced the lower limit of integration by because the integral is convergent for y — ► 0. Therefore, up to 
corrections suppressed by the inverse powers of cutoff A, Eq. (|B6|) is indeed equivalent to Eq. (|B5p . Then, by using 
the same approach as in the second paper in Ref. |28| . we expand the result on the right hand side of Eq. (|B6[) in 
powers of 1 /A and arrive at the following form of the gap equation: 
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where A = 4^4A/(y / 7re^) = G; nt A/ (Air z / 2 h 2 v 2 F ) is the dimensionlcss coupling constant and ((z,q) is the generalized 
Ricmann zeta function.— By assuming that the gap A s is much smaller than the Landau energy scale cb , we find the 
solution in an analytical form, 

A 3 =M=-^—, ji a = Aa±. (B9) 

1 — A 

Here, in order to get the result for the chemical potential p, s we used Eq. (|2"Tj) . It is easy to check that the gap equation 
also has another solution, which is obtained from Eq. (|B9[) by replacing A s and jl s with — A s and — jl s . However, the 
second solution is equivalent to that in Eq. (|B9|) : one can see this from dispersion relations (fl7|) . (| 18[) by transforming 
cr — > — a there. In other words, these solutions describe two degenerate ground states connected by a Z^s [C SU(2) S ] 
symmetry transformation. 

Turning to case (b) in Eq. (|B4[) . we have A s = 0, 

A s = ±M and p, s = 0, (BIO) 

where the last relation follows from Eq. (|2"Tj) . 

Finally, we would like to note that by analyzing the inequalities \/j, s £t s \ < \A S ± A s |, one can show that solution 
(|B9|) with a triplet Dirac mass exists for 

\H.\<M-A (Bll) 
and solution ()B10|) with a singlet Dirac mass exists for 

\H.\<M. (B12) 



2. The second case 

There arc two possibilities \^ s + fl s \ < |A S -A S |, \fx s -ft s \ > |A S +A S | or \[i a — pt a \ < |A S + A S |, \fi s +jl s \ > |A a — A a |. 
In the first case, the equations for Dirac masses take the form 

As + A a = -Aa ± sign(// fl - fl s ) + 2A A "^ As , (B13) 



n=l 



00 A _ A 

A s -A s = AJ2 V s [l + e(n-l)}, (B14) 

where s± = sgn(eB±). While the equation for A s — A s coincides with Eq. (|B2|) . the equation for A s + A s is slightly 
different from its counterpart in Eq. (|B1|) . Unlike Eq. (|B1[) . the above equation for A s + A s does not contain the sign 
factor sign(A s + A s ) in the LLL contribution. The absence of such a factor in Eq. (|B13|) means that the sign of the 
LLL contribution is fixed for a given set of values of ^t s , fi 8 , and eB±. In turn, this implies that Eq. (|B13[) [unlike the 
gap equation (|B1[) ] has only one solution whose sign is correlated with the sign of the LLL contribution. In order to 
prove this, let us consider the following equation: 

oo 

x = -A + 2AY ; f =■ (B15) 

By taking x negative, we see that its absolute value |a;| satisfies an equation that is equivalent to the equation for 
positive A s that follows from Eq. (|B5|) . Therefore, the solution for |x| coincides with the positive solution for A s in 
(|B9|) . We can also show that Eq. (|B15[) docs not have a solution for positive x by using the integral form of (|B15[) . 



a r d y 



-_e 



(B16) 



where the term —2 is subtracted in order to get the negative LLL contribution as in (|B15[) [cf. Eq. (|B6[) ]. 

In order to prove that Eq. (|B16[) does not have solution, we will use the fact that Eq. (|B6p does not have a 
nontrivial solution for B± — > in the case when the coupling constant G- ln t is subcritical, i.e., G; n t < 4n 3 / 2 Vph 2 / A, 
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or equivalcntly A < 1. Note that the coupling constant should indeed be subcritical because, as we know from 
experiment, there is no gap generation at B± = 0. It is not difficult to prove that the right hand side of Eq. (|B16|) 
is less than A after taking into account that i(cotht — 2) < I for t > 0. Then we conclude that Eq. (|B16|) does not 
have a solution for a subcritical coupling constant A < 1. As for Eq. (|B15|) . it has only one solution which, in fact, 
coincides with the solution for A s in (|B9|) times — 1. Thus, the solutions of Eqs. (|B13[) and (|BI4[) are 

A s + A s = -sign(/i s - ft s ) s ± M (BI7) 

and 

A s -A s =±(A S + A S ). (BI8) 

From the fact that the solutions for A s + A s and A s — A s have the same absolute value, we conclude that either 
A s 0, A s = or vice versa A s ^ 0, A s = depending on the sign in Eq. (|B18[> . If Eq. (|B16p had solution, there 
would exist solutions with both nonzero A s and A s . 

Further, solution of Eq. (|2ip for fl a in the case under consideration takes the form 

A 

fi s = — [-sign(/i s - fi s ) + sign(A s - A s ) s±\. (B19) 

Using (|B17p . it is easy to check that for the plus sign in (|B18[) (when A s ^ and A s = 0) fig = Asign(A s ) s± 
and for the sign minus in (|B18p (when A s = and A s ^ 0) fi s = 0. In the latter case, the assumed inequalities 
\fi a + fi s \ < \A S — A s \ and \fi a — fi a \ > \A S + A s \ cannot be satisfied, therefore, only solution with triplet Dirac mass 
A s is realized 

A s = — sx sign(/i s — fi s ) M, A s = 0, fi s = Asign(A s ) s±. (B20) 
In the other case \fj, a + jj, s \ > \A S — A s | and \fi s — fi s \ < \A S + A s |, we find the following solution: 

A s = s_l sign(/i s + fi s ) M, A s =0, ji s = Asign(A s ) s ± . (B21) 
One can show that it is possible to join solutions (|B20|) and (|B21|) into one solution with triplet Dirac mass 

A s = M, fig = As±, A s = 0, (B22) 

which exists for 

M — A< \fj, s \ < M + A. (B23) 

In fact, like in the previous subsection, there is another solution, with A s , fi s replaced by — A s , — fi s . However, such 
a solution is equivalent to solution (|B22[) by a SU{2) 8 (or Zi s ) symmetry transformation. 



3. The third case 

For ± fi s \ > \A S =F A s |, the equations for Dirac masses take the form 

A s + A s = -Asignfc, - fig) S± + 2AJ2 As ^ + As , (B24) 

n=l ^ ns 

°° A - A 

As-As = Asign(fi s +fis)s± + 2AY^ " (B25) 

n=l Lns 

Solutions of Eqs. (|B24)l and (|B25|) are 

A s + A s = -s_Lsign(/x s -fi s )M, (B26) 
A s - A s = s± sign(/i s + fig) M. (B27) 

Using these solutions and taking into account the inequalities \fi s ±/I s | > |A S =F A s |, one can check that Eq. (|2ip has 
only the trivial solution. Then it follows from Eqs. (|B26p and (|B27|) that 

As = -s± sign(^ s ) M, As = fis = 0. (B28) 

Taking into account the assumed inequalities |// s ±/2 s | > |A S =p A s |, we find that this solution with singlet Dirac mass 
exists for 

\fi s \>M. (B29) 
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4. Final solutions for A s , A s , and p, s as functions of fi a 

Using the results derived above and calculating the quantity X s in Eq. (|25[) (which is needed for solving the equation 
for fi s ), we obtain the following three distinct solutions. 

• Solution I (triplet Dirac mass). By joining the two solutions of the same type in Eqs. (|B9|) and (|B22j) . considered 
in Subsecs. IB H and IB 21 respectively, we arrive at the following solution: 

A S =M, Vs=As ± , A s = 0, X s = (B30) 

which exists over the combined range of validity \fx s \ < M + A. Let us mention that there is also another 
solution, in which A s and fl s are replaced by — A s and — fl s , respectively However, this second solution is 
equivalent to that in Eq. (|B30|) : one can see this from dispersion relations (fT"7j) . (fT8|) by transforming a — > —a 
there. In other words, the two solutions are related to two degenerate ground states connected by a SU(2) S (or 
i?2s) flavor transformation. 

• Solution II (singlet Dirac mass). This is one of the two solutions in Eq. (|B10[) from Subsec. IB II that corresponds 
to a particular choice of the sign for the singlet Dirac mass, 

A s = s x sign(/i s ) M, A s = fi s = 0, X s = AA sign(^ s ). (B31) 

It exists for \fi s \ < M. 

• Solution III (singlet Dirac mass). This combines the remaining solution in Eq. (|B10|) from Subsec. IB II with 
solution (|B28)) in Subsec. [B3] to give 

A s = -s± sign(^ s ) M, A s = fl s = 0, X s = -4Asign(/x fl ). (B32) 

This solution exists for all values of /i s . 

A noticeable point is that unlike the case with a triplet Dirac mass, the solutions II and III, with a different sign for a 
singlet Dirac mass, are different. This in particular can be seen from dispersion relation (fl8]l . This feature is directly 
connected with the fact that while the triplet mass is even under time reversal T, the singlet mass is T-odd. The 
latter is in turn connected with the fact that A oc s± = sign(B^) (recall that a magnetic field is also T-odd). 

Let us also emphasize that the expressions for Dirac masses in solutions I, II, and III are valid only for A < I: in 
the supercritical regime, with A > 1, a Dirac mass A is generated even with no magnetic fields Experiments clearly 
show that the subcritical regime, with A < 1, takes place in graphenc^ As argued in Sec. IIVI in the main text, 
realistic values for A in this model are A < 0.2. 



5. Including both spin up and spin down states 

In the previous subsection, the solutions for masses and chemical potentials were found for a fixed spin, treating 
the electron chemical potential fi s as a free parameter. Here we will describe full solutions, including both spin up 
and spin down states. For this purpose, we need to solve Eq. (|22[) for the chemical potentials [i±. Since the X term in 
that equation contains both spin up and spin down contributions, the equations for ji + and /i_ are now coupled and 
have to be solved together. As a result, the full chemical potentials [i± will be expressed through the bare electron 
chemical potentials p,± = /io T Z. 

At a fixed spin, there are 3 different types of solutions for masses and /2 S described in Subsec. IB 41 Since we can 
choose any of them for each spin, there are nine possible types and, therefore, nine systems of coupled equations for 
H + and Fortunately, noting that the solutions for the types II-I, III-I, and III-II can be obtained from those for 
I-II, I-III, and II-III by interchanging the spin subscripts + and — in the latter, this number can be reduced to six 
coupled systems. We will analyze them below case by case. 

It will be convenient to separate these systems of equations into 3 groups. The first group includes one system, I-I. 
This is the simplest case with triplet masses A± for both spins, when the Hartree diagram does not contribute in the 
equations for /j,± . The second group consists of hybrid systems I-II and I-III, where while fields with spin up have a 
triplet mass A + , the fields with spin down have a singlet mass A_. The third group, II-II, II-III, and III-III, consists 
of solutions with singlet masses A± only. 

In the analysis, it will be assumed that the Zeeman energy Z < A. As argued in Sec. IIVI this assumption is valid 
for magnetic fields \B±\ < 45T used in experiments! 13 ' 14 
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• The first group: Triplet Dirac masses 

o I-I. In this simplest case, using Eq. (|B30|) . we immediately find from Eq. (|22[) that /x+ = /2+ and the solution 
is: 

A+ = M, p, + = As x , M+ = P>+, A+ = 0, 
A_ = M, fi^=As_i_, [i-=p,-, A_ = 0. 

It exists for 

\p,+ \<A + M, \p,-\<A + M. (B34) 

The physical meaning of these constraints is clear: they imply that the LLL is neither completely filled nor 
empty. 

• The second group: Hybrid solutions 

o I-II. By using Eqs. (|B30[) and (|B31|) . we analyze the system of two equations |22p for /i+ and /«_ and find 
that the solution 



(B35) 



A+ = M, /2+ — As± , /x+ = /t+ — 4Asign(//+), A+ = 0, 
A_ = /2_ = 0, fX- = jl- — 3Asign(/2_), A_ = — s± sign(/2_) M 

exists for 

3A-M < 1/2+1 <5A + M, 3A - M < |/2-| < 3A, sign(/2+) sign(/2_) > 0. (B36) 
o I-III. In this case, using Eqs. l[2"2"]). (|B30|) . and (|B32|) . we find the solution 

A+=M, p+=Asi, M+ =M+ -4Asign(/Z + ), A+ = 0, 

(B3/j 

A_ = /}_ = 0, /t_ = /x_ — 3Asign(/i_), A_ = — sj_ sign(/t_) M, 

which exists for 

3A - M < 1/2+1 < + M, |/2_| > 3A, sign(/2+)sign(/2_) > 0. (B38) 

• The third group: Singlet Dirac masses 

o II-II. Using Eq. (|B31[) and analyzing equations ([22]) for /t+ and /t_, we find the solution 

A+ = /i+ = 0, /i+ = /!+ - 7Asign(/2+), A+ = -s± sign(/2+)M, g 
A_ = p,- =0, /j,- = fi- — 7 A sign(/2_ ) , A_ = — sj_sign(/2_)M, 

which exists for 

— M < 1/2+1 < 7A, 7 A - M < |/2_| < 7A, sign(/2+)sign(/2_) > 0. (B40) 
[Formally, there is also another solution, 

A+ = /2+ = 0, (i + = /2+ - Asign(/2+), A+ = s+ sign(/2+) M, 

(B41) 

A_ = /2_ =0, |U_ = /t_ — ^4sign(/i_), A_ = sx sign(/t_) M, 

which exists for 

A < 1/2+1 < A + M, A < |/2_| < A + M, sign(/2+)sign(/2_) < 0. (B42) 

However, because of the latter inequalities, it is easy to check that this solution does not satisfy the 
condition Z < A and therefore is not realized for magnetic fields \B±\ < 45T.] 
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o II-III. As in the previous case, there are two solutions. The first solution, II-III-l, 

A+ = /2+ =0, fx+ = /2+ - Asign(/2+), A+ = s± sign(/2+) M, 
A_ = /}_ = 0, /i_ = /2_ + A sign(/2 + ) , A_ = -s± sign(/2+) M 

exists for 



The second solution, III-III-2, is 



A+ = /i+ = 0, fx + = fx + + A, A + = -s± M, 
A_ = /2_ = 0, /t_ = /2_ - A, A_ = s_l M. 



It is realized for 



The third solution. III-III-3. 



A+ = /i+ = 0, /t+ = /t+ - A, A + = s_l M, 
A_ = /i_ =0, /t_ = /2_ + A, A_ = -s ± M 



takes place for 



(B43) 



A < 1/2+1 < A + M, /2_ sign(/2+) > -A. (B44) 

The second solution, II-III-2, 

A+ = /2+ = 0, fi + = /2+ - 7Asign(/2+), A+ = -s+ sign(/2+) M, 

(B45) 

A_ = //_ = 0, /«_ = /t_ — 7^4sign(/t_), A_ = — sj_ sign(/x_) M 

exists for 

7 A - M < 1/2+1 < 7 A, /2_ sign(/2+) > 7A . (B46) 

o III-III. There are three solutions in this case. The first solution, III-III-l, 

A+ = /2+ = 0, fx + =p + - 7Asign(/2+), A+ = -s x sign(/2+) M, 
A_ = p,- = 0, /t_ = /2_ — 7Asign(/2_), A_ = — s± sign(/2_) M 

exists for 

1/2+1 >7A, |/2_|>7A, sign(/2+)sign(/2_) > 0. (B48) 



(B49) 



ft + >-A, /2_<A. (B50) 



(B51) 



p+<A, p_>-A, (B52) 
i.e., in fact, it is obtained from the second solution by interchanging spin subscripts + and — . 

6. Dependence of solutions on electron chemical potential po and free energy density of their ground states 

The process of filling LLs is described by varying the electron chemical potential /iq. Therefore, it will be convenient 
to express the intervals of the existence of the solutions found in the previous subsection in terms of /to- Henceforth 
we will consider no > 0. (Dynamics with negative /io is related by electron- hole symmetry and will not be discussed 
separately.) 

Some intermediate results of the analysis in this subsection will depend on whether the inequality M > 2Z or 
M < 2Z is satisfied. We will consider both these cases and indicate explicitly which inequality is satisfied for 
a particular solution. If nothing will be said, this means that the corresponding results are valid in both cases. 
Fortunately, the final results do not depend on whether M > 2Z or M < 2Z. 
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TABLE I: Intervals of the existence of solutions, relevant for the dynamics in the LLL at T — 0. 





A/ > 2Z 






M < 2Z 


T T 
1-1 


0<Ho<M + A-Z 






u \ Ho ^ m ~r zi — zj 


T TT 

i-ii 


M + A+Z<ho<3A- 


- Z 




no solution 


TT T 

11-1 


3 A - M + Z < no < 3 A + Z 




o/i — in 4- zv <. ho \ *jJ\ i zj 


T TTT 
1-111 


3 A - Z < ho < 5 A + M + Z 




1 A Aft X- 7 S ^ ^ A X. A T 4- 7 


HI T 
111-1 


3 A + Z < ho < 5 A + M 


- z 






TT TT 
11-11 


7 A - M + Z < no < 7 A 


- z 




no solution 


TT TTT 1 
11-111- i 


A + Z < no < M + A + Z 




J± -f- Zj <, /1(J <v iVi 4- /I 4- Zj 


TT TTT 
11-111-Z 


7 A 7 nr. 7 A 4- 7 

l -rl — Zj \ [AO ~~ i/1 t ^ 






7 A AyT 1 7 ^ . . . ^ ^7 /I 1 7 


iii-ii-i 


A 7 iir, <^ AT 4-4 

/l Zj \ fJ>0 i" T ^ 


z 




A-Z</i <M-|-yl-zT 


III-II-2 








no solution 


III-III-l 


74 4- 7 ^ //n 






7^4 + Z < ho 


III-III-2 


U < . /it) ' Zj 






< ho < A - Z 


III-III-3 


f) <T //n <^ 4 4- 7 






< ho < A + Z 


TA DT L 1 TT 

at T = 0. 


: The list of solutions that coexist in a set of non 
Ihe solutions with the lowest tree energy density 


-overlapping intervals of ho, relevant for the dynamics in the LLL 
are marked by stars. 


# 

77 


Interval 




M > 2Z 


M < 2Z 


1 


0<Ho< A- Z 




I-I, III-III-2, III-III-3* 


I-I, III-III-2, III-III-3* 


o 
z 


A-Z<ho<A + Z 




I-I, III-II, III-III-3* 


I-I, III-II, III-III-3* 


Q 
O 


A + Z < ho < M + A — Z 




I-I, III-II, II-III-l* 


I-I, III-II, II-III-l* 


A 
-± 


M + A-Z<ho<3A-M + Z 




II-III-l* 


II-III-l* 





3 A - M + Z < ho < 2A + Z 




I-II, II I, II-III-l* 


I-III, II-I, II-III-l* 


U 


2A + Z<ho<M + A + Z 




I-II*, II-I, II-III-l 


I-III*, II-I, II-III-l 


7 
t 


M + A + Z<ho<3A~Z 




I-II*, II-I 


I-III*, II-I 


o 
o 


3 A - Z < ho < 3 A + Z 




I-III*, II-I 


I-III*, II-I 


Q 


3A + Z < ho < 5A + M — Z 




i-iir, ni-i 


I-III*, III-I 


J.U 


5 A + M - Z < ho < 7 A - M + Z 




i-iii* 


I-III* 


ii 


7 A - M + Z < ho < GA + Z 




i-iii*, ii-ii 


I-III*, II-III-2 


12 


6A + Z < ho <5A + M + Z 




i-m, ii-ii* 


I-III, II-III-2* 


13 


5 A + M + Z < ho < 7 A - Z 




ii-ii* 


II-III-2* 


14 


7 A - Z < ho < 7 A + Z 




II-III-2* 


II-III-2* 


15 


7 A + Z < ho 




iii-iii-i* 


III-III-l* 



Taking into account that jl± = fio =F Z, we find the intervals of existence for solutions. These are given in Table [U 
Using this information, we see that some solutions may coexist. The list of coexisting solutions for a set of non- 
overlapping intervals of fio is summarized in Table UTl [We assume that Z > M — A = A\/(\ — A) which is likely to 
be satisfied because, as will be shown in Sec. IIV1 realistic values for A in this model are relatively small, A < O.2.] 

Thus, there are several coexistent solutions on different intervals of Ho- hi order to find the most stable solution 
among them, we have to calculate the free energy density Q of the ground states corresponding to these solutions. 
To facilitate this, we first calculate the free energy densities of the fixed spin solutions I, II and III considered in 
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Subscc. IB 41 by using expression (|C19[) for f2 derived in Appendix [Cl The results are 

solution I: fii = - !^±I [M + A + h], (B53) 

Anne 

solution II: fi n = - i^^i [A/ - (p + /2)sign(/i) + /il , (B54) 

solution III: » m = - [M + (/i + /2)sign(/i) + /t] , (B55) 



where h is the higher LLs contribution, defined by 



2M 4 M q 
~ 2d 



«=i ^/nel + A/" 2 fV^I + 71/2 + V^ £ -b) 



. 3\ , /5\ M 2 „ / M 4 



(B56) 



where is the Riemann zeta function. On the right hand side we used the expansion in powers of (M/cb) 2 - When 
keeping only the first two terms in the expansion, we find that the result deviates by less than 1% from the exact one 
for M < OAeg. Note that the above contribution from higher LLs is the same for all solutions. Therefore, it is only 
the LLL contribution that is relevant for choosing the lowest free energy density. 

It is not difficult now to calculate the free energy densities for all the solutions. In Table UH the solutions that have 
the lowest values of f2 and thus correspond to the ground states in the given intervals of /zo are marked by stars. As 
for the explicit expression for the energy density in the ground state, it reads 

n = _l£^M [M + A + 2Z+h], for 0</j, <2A + Z, (B57) 
zirnc 

V = --^^-[M -A + Z + h + Ho], for 2A + Z < ^ < + Z, (B58) 
zirnc 

^ = _l£#xl \m -7A + h + 2no], for 6A + Z<pt , (B59) 



Using now the explicit form of the solutions obtained in Subscc. IB 5[ wc can significantly reduce the number of the 
cases. As result, wc conclude that only the following three solutions are realized: 

(i) The solution with singlet Dirac masses for both spin up and spin down: 

A+ = LL+ = 0, LL+ = LL+ - A, A + = S ± M, 

. (B60) 
A_ = //_ = 0, fi- = /!_ + A, A_ = -s±M. 

It is the most favorable for < < 2A + Z^ We will call it the SI solution, which is one of several solutions 
with nonvanishing singlet Dirac masses. 

(ii) The hybrid solution with a triplet Dirac mass for spin up and a singlet mass for spin down: 

A+ = M, fi + =As ± , m+=M+-4A, A+ = 0, 
A_ = /}_ = 0, (i- = fi- — 3A, A_ = -sj_M. 

It is the most favorable for 2A + Z < [Iq < 6A + Z. Wc will call it the HI solution. 

(iii) The solution with equal singlet masses for both spin up and spin down: 

A+ = /i + = 0, n+ = /!+ - 7 A, A + = -s±M, 
A_ = /i_ = 0, fx- = fi- - 7A, A_ = -s±M. 

It is the most favorable for /j.q > 6 A + Z. We will call it the S2 solution. 



(B61) 



(B62) 
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APPENDIX C: FREE ENERGY DENSITY 



In this Appendix, the units with h = 1 and c = 1 are used. 

In order to calculate a free energy density f2, it is convenient to use the Baym-Kadanoff formalism (the effective 
action formalism for composite operators) developed in Ref. [55| (see in particular the last paper there). In the mean 
field approximation that we use, the corresponding effective action T has the following form: 

r(G) = -iTr[LnG- 1 +5- 1 G-l] +^ J d 3 x {tr [-y°G(x, x)7°G(a;, x)] - (tr [-y°G(x, x)\ f } , (CI) 

where the trace, the logarithm, and the product S~ 1 G are taken in the functional sense, and G = diag(G+, G_). The 
free energy density f2 is expressed through T as Q = — T /TV, where TV is a space-time volume. The stationarity 
condition ST(G)/5G = leads to the gap equation (JSJ. On its solutions we have 



-i Tr 



LnG -1 + - [S- X G- 1) 



(C2) 



Henceforth we will use the symmetric gauge with A(r) = (— Bj_y/2, Bj_x/2). Then, as was shown in Appendix [Al 
the Green's function G s (u, u'), with u = (t, r), has the form: 



G s (u,u')=e^ u ' u '^G s (u-u>), 

where $(w, u') = — er ■ A(r') is the Schwinger phase in the symmetric gauge. 
Because of the translation invariance in time, we have 



(C3) 



f ^e-^-*')G s (c;r,r'). 



(C4) 



Then the effective action T can be rewritten as 



, dui 
-iT / — Tr 
2tt 



lnG" 1 ^) + - (5" 1 (cj)G(cj) - 1) 



(C5) 



where 



G a 1 (u;r i i f ) = [(cj + A i s ) 7 ° -^(tt -7) +i/i s7 1 7 2 +«A S 7°7 1 7 2 
^- 1 (^;r,r') - -i [(w + }l s ) 7 - v F (n ■ 7 )] S(r - r'). 



S(r-r'), 



(C6) 
(C7) 



In Eq. (|C5[) the functional operation Tr includes now the integration over the space coordinates only and the trace 
over matrix indices. 

Integrating by parts the logarithm term in Eq. (|C5[) and omitting the irrelevant surface term (independent of the 
physical parameters), we arrive at the expression 



-iT 



duj r 
2^ 



-Tr 



-w^^-G(w) + ±(S-\u,)G(u)-l) 



(C8) 



with 



OG-^lo) 



-ij°6(r-r'). 



(C9) 



Substituting now expression (|C3[) for the Green's function in T, one can see that the Schwinger phase goes away and 
we get 



T = -iTV 



dw 
2tt 



-tr 



«7°wG( W ; 0) + 1 (-i [(w + /2)7° - ■ 7 )] G{u; r)| r=0 - 5(0)) 



(CIO) 
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Dividing T by the space-time volume TV, we find the free energy density: 







duj f d 2 k 



2ir J (2tt) 2 
dw f d 2 k 



4tt J (2tt) : 



tr \ iuyfG(w, k) + - (-i [(to + ^) 7 ° - t- F (k • 7)] k) - l) 



■tr { [(w - m)7° + ■ 7)] k) + 



(Cll) 



where the propagator G s (u,k) is given in Eq. (|A21jl in Appendix lAl By making use of its explicit form, we can 
calculate the following two integrals that contribute to the free energy density, 



d 2 k 



7°e.(*,k) = t^E 



(u + Us + ij^sl ! 1 ! 2 - *A a 7 1 7 2 + A s7 °) P„ 



ri 2 fc 
(2tt)< 



« F (k. 7 )G s Kk) = -Lg 



v 2 F \eB ± \n9(n-l) 



nP t'o ^ + /*- + ^7°7 X 7 2 ) 2 - (A s - iA fl7 V7 2 ) 2 - 2v 2 F \eB ± \i 



where 



P n = 1 - i 7 1 7 2 sign(eB_ L ) + [l + i 7 1 7 2 sign(eS ± )] 0(n - 1) 
In the calculation, we used formula 7.414.7 from Ref. (EH, i.e., 



e- at t a Ll{t)dt 



T(a + n + l)(a-l) n 



n'.a 



\ n a+n+l 



, Re a > — 1, Re a > 0. 



(C12) 



(C13) 



(C14) 



(C15) 



By dropping an infinite divergent term which is independent of the physical parameters, from Eq. (|C11[) we derive 
the following expression for the free energy density: 



o = - 



(4^) 2 *Z 



B = ± 



E / dwtr ^E 



(UJ ~ fig) 



uj + fi s + ip.s'f^Y ~ i&sTT + A si 



P n +Av 2 F \eB ± \n9(n- 1) 



(w + lis + iHs-y ! 1 ! 2 ) 2 - (A, - zA s7 o 7 i 7 2 ) 2 - 2v 2 F \eB ± \ 



(C16) 



Here the trace trp is taken over the Dirac indices. 

The free energy density O is a function of A s , fx s , /i s , A s , p, s , and Bj_. Normalizing by subtracting its value at 
A s = pj S = (i s = A s = £L S = 0, we obtain: 



v ' s=±n=0 *L 



(h> - p. a )[u + Ms + tAs7°7 1 7 2 ~ *A fl 7V + A s7 °]P„ + Av 2 F \eB ± \n9(n - 1) 
(cj + iesign(w) + (j, s + ^ s7 7 1 7 2 ) 2 - (A s - zA s7 ° 7 1 7 2 ) 2 - 2v 2 F \eB±\n 



Lu 2 P n +4v 2 F \eB ± \n9(n - 1) 
(uj + iesign(cj)) — 2v F \eB±\n 



One can check that for jl s = A s = fi s = fi s = B± = and A + = A_ = A this expression reduces to 



(C17) 



O(A,0,0,0,0,0) = -^ 



dx 



_A3 
2 ~ 6tt : 



v / A 2 + x(v / A 2 + .t + % /^) 
which coincides with the known expression for the vacuum energy density in 2 + 1 dimension 



(C18) 
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Finally, integrating over lo and taking trace, we find the following expression for the free energy density: 



= _ g^2 \ [a*s + Ms - Ms - (A a + A s )sign(eB ± ) sign(/i s - p, s )9{\^ s - p,„\ - |A S + A s 
+ [A s + A, - + Ms - /2 s )sign(eB ± ) sign(A s + A s )0(|A s + A s | - \fi a - /i s |) 



+2E 



+ Ms - Ms)sign(// S - - 2e B \pn\ (|m s - - 



(A, + A s ) 4 0(ff+ - | Ms - 



[Ms -> -Ms, A s -» -A s , sign(eB ± ) -» -sigii(eBj_)] 



(C19) 



where ££f. 



(A s ± A s ) 2 and e B = v^leBj 



APPENDIX D: ANALYTIC SOLUTIONS OF GAP EQUATION FOR n = 1 LL AT T = 

1. Fixed spin 



In Appendix [B] we analyzed solutions of the gap equations under the condition that only states on the LLL can be 
filled, |ms ± /is | <C €b = y/^\sB±\vp/c. Since all the dynamically generated parameters are much less than e B , this 
condition implies that the bare chemical potential mo also has to satisfy mo "C £_b in that case. 

In this section, we will consider the case when mo is of the order of the Landau scale e B , i.e., we will study the 
dynamics when states on the first Landau level, n = 1 LL, can be filled. The gap equations are given in Eqs. (|19|) -([22 |) 
in Sec. IIIII In order to get their solutions for mo ~ € b, wc will follow the steps in the analysis in Appendix IBl The 
equations for the dynamical parameters A s , A s , and jl s form independent systems of equations for each spin. From 
these systems, we can find their solutions as functions of Ms- Wc obtain the following three solutions. 



• Solution f-I. This solution corresponds to the case with |m« — Ms I < + (A s + A s ) 2 and |m« + M«| < 



(A, - A s ) 2 . It is 



A s = -sj_ sign(M s ) M, A s = Ms = 0. 



(Dl) 



This solution exists for M < \/j, s \ < y/e% + M 2 ■ Actually, it is exactly the same as solution (|B28j) considered in 
Subsec. 3 of Appendix [Bl For positive Ms, this solution corresponds to a state with the completely filled LLL 
and the empty n — 1 LL. With increasing Ms, this solution exists up to the point where the first LL starts to 
fill, which is defined by the upper limit of the above inequality for Ms- Recall that X s = — 4Asign(M s ) for such 
a solution. 



Solution f-II. This solution is realized when the inequalities \fi a — jl s \ < \J e 2 B + (A s + A s ) 2 and |m 



> 



e 2 B + (A s — A s ) 2 are satisfied. In this solution, all three dynamical parameters A s , A s , and fi s are nonzero: 



-sxsign(M s ) 



AI - Ah 



A, 



■s± sign(M s ) 



M 1 + AI 



M s = sign(M s ) A. 



(D2) 



Here Mi satisfies the following equation: 



1 = 



.4 



/A? y/V 



-yMl 



cothj ^-y) -2exp(-ye 2 B ) 



(D3) 



Note that the last term in the square brackets of the integrand appears because the n = 1 LL contribution is 
absent in the equation for A — A [cf. Eq. (|B6|) where all LLs are included] . 

Utilizing the analysis in the second paper in Ref. [28l . we arrive at the following gap equation for M\ : 



A 2A/1 M? 
Mi e B \2 ei 



2A 



V4 + 



O A 



A 2 " 



(D4) 
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where £(z, q) is the generalized Riemann zeta function^ In the subcritical regime (A < 1) its solution is given 

by 

M > * l-A + 2[l-<(l/2)]A/ eB ' (D5) 

Since the last term in the denominator is positive, we have M\ < M that is consistent with the fact that the 
equation for A — A misses the contribution of the n = 1 LL. 



This solution exists for \/e% + M\—A < \fj, s \ < + M 2 +A. One can check that the corresponding parameter 
X s is X s = — 8Asign(/i s ). As in the case of the LLL (see Subsec. IB 4p . there is another solution with A s , fb s 



replaced by — A s , —jj, s , which takes place for \fj, s —p, s \ > y e 2 B + (A s + A s ) 2 and |/i s +/2 S | < y e B + (A s — A s ) 2 . 
These two solutions are equivalent: one can see this from dispersion relations (|17[) . (fT5|) by transforming a — » — a 
there, i.e., as in the case of the LLL solution I (|B30[) . these solutions are related to two degenerate ground states 
connected by a SU(2) S (or Z^s) flavor transformation. 



• Solution f-III. This solution corresponds to the case with \fi s — /2 S | > y e 2 B + (A s + A s ) 2 and \fi s + fl s \ > 
e 2 B + (A s — A s ) 2 . Its explicit form reads 

A s = Jx s = 0, A s = -s± sign(/x s )Mi. (D6) 



This solution takes place for |/z s | > \J t B + M( and the corresponding X s is X s = — VIA sign(/x s ). 



2. Including both spin up and spin down 

In Subsec. TP 1( the solutions for the dynamical parameters A s , A s , and jl s at fixed spin were described. Since X 
contains contribution of fields of both spins, the equations for chemical potentials fj,+ and /j,- for fields of different 
spin are coupled and have to be solved together. Since we can choose any of the found three solutions for masses at 
a fixed spin, we should solve 9 systems of coupled equations for /i + and Like in the case of the LLL, it is enough 
to consider only 6 systems. The simplest case is the solution f-I-f-I because it corresponds to the case of completely 
filled LLL, which was already considered in Subsec. IB 51 We have 



f-I-f-I solution is given by 



A+ = /i+ = 0, /i+ = /!+ - 7Asiga(n + ), A+ = -s± sign(^ + )M, 
A_ = /i_ =0, /i_ = fi- — 7^4sign(/2_), A_ = — s± sign(/2_)M. 



(D7) 



It exists when sign(^ + ) sign(/j,_) > and 



7A + M < |/2 + | < 7A+ y e | +M 2 , 7 A + M < |/2_| < 7 A + y 1 e B + M 2 . (D8) 

It coincides with the solution III-III-l in Eq. (|B47|) in Subsec. IB 51 except for having a different lower limit for 
|/z±|. The latter is connected with the point that while the solution III in Eq. (|B32j) exists for all values of fi s , 
the solution f-I in Eq. (|D1|) exists only for |^ s | > M. This is because, according to the analysis in Subsec. IB 4| 
the solution III is a combination of two solutions: the solution (|B28[) . which is equivalent to the solution f-I, 
and one of the two solutions in Eq. (|B10|) . 

f-I-f-II solution is given by 

A+ = /i+ = 0, fx+ = /!+ - llAsign(/2 + ), A + = -s± sign(/2+)M, 

~ M-M x A _ t . ,_ . . ,_ ,M + Mi (°9) 
A_ = , [i- = -Asj_, = /!_ - 10^sign(/i_), A_ = -sj_sign(^_) . 

It exists when sign(/2 + ) sign(/2_) > and 

11A + M < 1/2+1 < UA+ y/e 2 B +M 2 , 9A + y^^M 2 < \fi_ | < 1L4 + yU^+M 2 . (D10) 
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f-I-f-III solution reads 



A+ = /2+ = 0, = - 15Asign(/i + ), A + = -s ± sign(/i + )M, 

A_ = /i_ = 0, /i_ = p,- — 13^4sign(/I_), A_ = — s± sign(/2_)Mi, 



and takes place when sign(/x_|_) sign(/z_) > and 



(Dll) 



15A + M < \p, + \ < 15 A + y + M 2 , 13 A + ^ e% + M\ < \p._ |. (D12) 

f-II-f-II solution is given by 

~ M - M\ . _ . . » . , .M + Mi 

A+ = , fj,+ = -As±, n+ = n+ - UAsign(jj, + ), A+ = -s± sign(/i+) , 

s M-Mi _ , . . ._ ,M + Mi ^° 13 ^ 

A- = , /i- = -As_l, /i_ = /x_ - TL4sign(jU_), A_ = -s± sign(/z_) , 

and exists when sign(/2 + ) sign(/2_) > and 



13A+ Je| +M? < |/2±| < 15A + t/4 + AP - ( D14 ) 



f-II-f-III solution is given by 



- M -Mi „ _ ._ . . , .M + Mi 

A+ = , /U + = -Asx, ^+ = ~ 18Asign(/i + ), A+ = -s_LSign(^ + ) , 

A_ = /i_ =0, jU_ = /2_ — 17Asign(/2_), A_ = — sj_ sign(/i_)Mi, 
and takes place when sign(/2 + ) sign(/i_) > and 



17A+Je%+M? < l/i+l < 19A+ Je%+M 2 , \fi_\ > 17 A + J e B + M\ . (D16) 



• f-III f-III solution is given by 



A+ = p.+ = 0, ^+ = n+ - 21Asign(/i+), A + = -s\_ sign(^+)Mi, 
A_ = /i_ = 0, /i_ = //_ — 21^4sign(/2_), A_ = — sj_ sign(/z_)Mi, 



and exists when sign(/i + ) sign(^_) > and 



(D17) 



\ft + \>21A+y l e B +M(, \p.\>21A + y/e B +Mf. (D18) 

3. Dependence of solutions on [Mo and their free energy density energy 

Using the solutions found in the previous subsection, we find that the intervals of their existence in terms of /iq for 
A*o > (dynamics with negative /io is related by the electron- hole symmetry and will not be discussed separately). 
These are given in Table IIII1 By making use of this information, we can also determine the complete set of non- 
overlapping intervals of fio and the solutions that (co-)exist on such intervals. This is summarized in Table ITVl 

Thus, there are several coexistent solutions on certain intervals of fiQ. In order to define which solutions arc realized, 
we have to calculate their free energy densities. To facilitate this calculation, first we will calculate free energy densities 
of solutions f-I, f-II, and f-III. Using the effective potential given by Eq. (|C19[) . we have 

solution f-I: n H = - [M + (a + a) sign(u) + h] , (D19) 

4irnc 

solution f-II: Q UI = ^£±Ml -\- A + 2(/u + ft) sign(/x) - 2e B + , (D20) 

solution f-III: fi f _ m = {M x + 3(fx + ft) sign(/i) - 4e B + hi) , (D21) 
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TABLE III: Intervals of the existence of solutions, relevant for the dynamics in the n = 1 LL at T = 0. 



f T-f T 
i-i i-i 


7 A 4- 4- 7 f //,. ^ 74 x . /f 2 4- M2 7 

I Z± ^ 1V1 T Zj flQ < 1 zi -w fc£j ill zj 


f T-f TT 
i-i i-ii 




f TT-f T 
i-ii i-i 


Q 4 4- . 1 e 2 4- A f 2 4- 7 <r- lir. <r 1 1 4 4- . //^ 2 4- A f 2 7 
yzi -(- y -}- -(- Z/ ^ /io - 11/1 -f- fc£j 4^ ill Zj 


f-I-f-III 


H/lx, / f 2 4_ A,/" 2 7 <- ,,„ ISi X, / f 2 _L Af 2 4- 7 


f-III-f-I 


13,4 + ^4 + A/f + Z < /xo < 15A + ,/e| + Af 2 - Z 


f-II-f-II 


13,4 + + + Z < no < 15A + ^/e| + M 2 - Z 


f-II-f-III 


17A + + Mj 2 + Z < no < 19A + ^4 +M 2 + Z 


f-III-f-II 


17 A + + M'{ + z7 < no < 19A + y/e 2 B + M 2 - Z 



f-III-f-III no > 21A + + M 2 + Z 



TABLE IV: The list of solutions that coexist in a set of non-overlapping intervals of ^o, relevant for the dynamics in the n — 1 
LL at T = 0. The solutions with the lowest free energy density are marked by stars. 



# 


Interval 


Solution(s) 


1 


7 A + M + Z < no < 7A + ^/e% + M 2 — Z 


f-I-f-I* 


2 


9A + + M 2 - Z < no < 9A + v/e| + M 2 + Z 


f-I-f-II* 


3 


9A + y/e 2 B + M( + Z < no < HA + ^4 + M 2 — Z 


f-I-f-II*, f-II-f-I 


4 


13A + y/e B + Ml - Z < no < 13A + V4 + M i + z 


f-I-f-III* 


5 


13A + yjt\ + Ml + Z < no < ISA + s/e\ + M 2 - Z 


f-I-f-III*, f-III-f-I, f-II-f-II 


6 


15A + ^/e 2 B + M 2 - Z < no < 15A + y/e* B + M 2 + Z 


f-I-f-III* 


7 


17 A + y/ej, + Ml + Z < no < 19A + V e l + M 2 - Z 


f-II-f-III*, f-III-f-II 


8 


19 A + y/e B + M 2 - Z < no < 19 A + + M 2 + Z 


f-II-f-III* 


9 


no > 21A + y/e 2 B + Ml + Z 


f-III-f-III* 



where h is given in Eq. (|B56|) and 



2Mf 



r.- . 

n=2 y/ n e 2 B + Mf [yJn€ 2 B +Ml + y/ne B 



M 4 r 



2 ~ 2e 3 



M 2 / M 4 \ 
C(3/2)-l-[C(5/2)-l]-^- + 0(-^ 



(D22) 



Now it is not difficult to calculate free energy densities for all solutions and determine the ground state on each 
interval. The solutions with the lowest free energy density are marked by stars in Table HVl The explicit form of the 
corresponding energy densities are 



f-I-f-I: 
f-I-f-II: 

f-I-f-III: 

f-II-f-III: 

f-III-f-III: 



O 
9. 
O 

n 



2-kTlc 



eB, 



2nhc 
eB ± 



2nhc 
eB ± 



2nhc 
eB ± 



(M + 2/i - 7 A + h) , 
/ 3M + Mi 



V 4 
fM + Mi 

/ 2>M\ + M 



-3no-^A-e B + Z- 
Ano~27A-2e B + 2Z 
-5Mo-43A-3e B + Z 



2nhc 



(Mi + 6^o - 63A - 4e B + hi) . 





(D23) 


3h + hi \ 
4 J' 


(D24) 


h + hA 
2 J' 


(D25) 


3hi +h\ 
+ 4 )• 


(D26) 
(D27) 



Therefore, the number of different solutions is reduced down to following five. 
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(f-i) The solution f-I-f-I 

A+ = fi+ = 0, n+ = n+ - 7 A, A+ = -s± M, 



(D28) 

A_ = /i_ =0, fi-=p-- 7 A, A_ = -s_l M 



is realized for 7 A + M + Z < /i < 7 A + yJe 2 B + M 2 — Z and has free energy density in Eq. (|D23|) . This result 
means that the solution S2 given by Eq. (|B62[) in Subsec. IB 61 takes place for /io < 7 A + \J ( 2 B + M 2 — Z. 

(f-ii) The solution f-I-f-II 

A + = p,+ = 0, n + = p+ - 11 A, A+ = -s± M, 

x M-Mi . , . . M + Mj (D29) 

A_ = , /i_ = -Asj_, /Lt_ = /i_ - 10A, A_ = -s± 

takes place for 9 A + + M| - Z < fi Q < 11 A + i/ej + M 2 - Z and has free energy density in Eq. (|D24j) . 
(f-iii) The solution f-I-f-III 

A+ = ii+ = 0, u+ = fi+ — lbA , A + = — sj_ M, 

„ + + + (D30) 

A_ = /i_ = 0, = /Z_ - 13.4, A_ = -s_l Mi 



is realized for 13,4 + ^e 2 , + M 2 - Z < // < 15A + y/e 2 B + M 2 + Z and has free energy density in Eq. (|D25j) . 

(f-iv) The solution f-II-f-III 

~ M-Ah A _ M + Mi 

A+ = , ^+ = -As+, n+ = n+ - l&A, A+ = -sj_ - , ^ D3i ^ 

A_ = /i_ = 0, = /t_ - 17A, A_ = -s± Mx 



takes place for 17 A + y / e|+Mf + Z < /i < 19^4 + y/e B + M 2 + Z and has free energy density in Eq. (|D26|) . 
(f-v) The solution f-III f-III 

A+=fL+ = 0, M+ =/2+-2L4, A+ = -8±M 1 , 
A_ = /i_ = 0, /i_ = /i_ - 21A, A_ = -s ± M-y 



is realized for /iq > 21 A + y/e B + Mf + Z and has free energy density in Eq. (|D27I) . 
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